Stability of general-relativistic accretion disks 
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Self-gravitating relativistic disks around black holes can form as transient structures in a number 
of astrophysical scenarios such as binary neutron star and black hole-neutron star coalescences, as 
well as the core-collapse of massive stars. We explore the stability of such disks against runaway 
and non-axisymmetric instabilities using three-dimensional hydrodynamics simulations in full gen- 
eral relativity using the THOR code. We model the disk matter using the ideal fluid approximation 
with a T-law equation of state with T = 4/3. We explore three disk models around non-rotating black 
holes with disk-to-black hole mass ratios of 0.24, 0.17 and 0.11. Due to metric blending in our initial 
data, all of our initial models contain an initial axisymmetric perturbation which induces radial disk 
oscillations. Despite these oscillations, our models do not develop the runaway instability during 
the first several orbital periods. Instead, all of the models develop unstable non-axisymmetric modes 
on a dynamical timescale. We observe two distinct types of instabilities: the Papaloizou-Pringle and 
the so-called intermediate type instabilities. The development of the non-axisymmetric mode with 
azimuthal number m = 1 is accompanied by an outspiraling motion of the black hole, which signifi- 
cantly amplifies the growth rate of the m = 1 mode in some cases. Overall, our simulations show that 
the properties of the unstable non-axisymmetric modes in our disk models are qualitatively similar 
to those in Newtonian theory. 

PACS numbers: 



I. INTRODUCTION 



Thick relativistic accretion disks and tori around black 
holes (BHs) can form as transient structures in several 
astrophysical scenarios, including core-collapse of mas- 
sive stars [1-12 ] and the merger of neutron star (NS) 1 13- 
|2Zl and NS-BH binaries d |28H37||. Many models of 
gamma-ray bursts (GRBs) rely on the existence of mas- 
sive dense accretion disks around BHs [1.3, 38-40]. The 
observed ~ 10 51 erg energy powering GRBs [39 4T1I421 
is believed to be coming either from the accretion disk 
and/or rotation of the central object. If this energy 
comes from the disk 1 , then - assuming that the effi- 
ciency of converting disk energy into that of GRB can 
at most be ~ 10% as in many other astrophysical sce- 
narios [46] - the accretion disk should have a mass 
of > 0.01 Mq. Recent numerical simulations have 
demonstrated that the mass of the disk resulting from 
binary NS-NS (BH-NS) mergers, which are thought to 
be candidates for central engine of short GRBs (see, 
e.g. 11401 |47| , but also 114811 ), can be in the range of ~ 
0.01 - 0.2 Mq m ED 123 El- Due to their larger ini- 
tial mass, in core-collapse of massive stars, which are 



believed to be progenitors of long GRBs 2 QHHH, signif- 
icantly more mass is available for forming disks 13l l50l . 

The neutrino-annihilation mechanism for triggering 
GRBs [40 51 52 j, in which neutrinos emitted by the disk 
annihilate predominantly at the rotation axis to produce 
e + — e~ pairs and deposit energy behind the jet, rely on 
super-Eddington accretion rates, which can take place 
only in high-density (p ~ 10 12 g cm~ 3 ) disks. Moreover, 
the efficiency of neutrino annihilation at the rotation axis 
and the requirement of a small baryon load of the rela- 
tivistic ejecta [39 ,41 , 53 1 was shown to strongly favor a 
toroidal structure for accreting matter 1541 155] . 

The duration of prompt 7-ray emission is > 2 s (< 2 s) 
for long (short) GRBs, while that of the later-time non- 
7 ray emission can be as long as ~ 10 5 s (e.g., I42T ). 
Recent observations of GRB afterglows have revealed 
a variety of late-time emission processes, including X- 
ray plateaus, flares, and chromatic breaks [42], some of 
which can persist up to ~ 10 5 s following the initial 
GRB prompt emission (see, e.g., [42 ] for a recent review). 
The amount of energy released in the late-time emis- 
sion phase can be comparable or even larger than that 
produced during the prompt emission phase. Both the 
prompt and the late-time emissions can be explained as 



In alternative models, the GRB is powered by the rotation energy 
of the central object. If the latter is a BH, then the rotation can be 
converted through Blandf ord-Znajek mechanism ] 43 1 . If the central 
object is a NS, then its energy can be transformed through magnetic 
fields l44l[45) . 



2 Note that although all of the long GRBs are believed to be produced 
by core-collapse of massive stars, not all of the latter can produce 
long GRBs: In order to produce a GRB, precollapse star is probably 
required to be rapidly rotating 1 1 12 49]. 
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the result of the activity of the central engine (see, e.g., 
Il42l l56l ), although alternatives models exist (see, e.g., 
II57H59I ). If the emission energy comes from the disk, 
such long durations of observed GRB emission either re- 
quire the disks to accrete in a quasi-stable manner for a 
sufficiently long period of time, or require the engine to 
be restarted in some way. 

Early studies of the stability of accretion disks have 
revealed that they can be subject to several types of ax- 
isymmetric and/or non-axisymmetric instabilities in a 
number of circumstances [60—66 [. Instabilities can lead 
to highly variable and unstable accretion rates, posing a 
serious challenge to the viability of "accretion-powered" 
GRB models. The so-called dynamical runaway instabil- 
ity of thick accretion disks around BHs was first discov- 
ered by Abramowicz, Calvani & Nobili [60 J. This insta- 
bility is similar to the dynamical instability in close bi- 
nary systems, when the more massive companion over- 
flows its Roche lobe. In such a case, the radius of the 
Roche lobe shrinks faster than the radius of the com- 
panion, leading to a catastrophic disruption of the latter. 
In disk+BH systems, a toroidal surface analogous to the 
Roche lobe can be found. A meridional cut of this sur- 
face has a cusp located at the I4 Lagrange point. If the 
disk is overflowing this toroidal "Roche lobe", then the 
mass-transfer through the cusp will advance the cusp 
outwards inside the disk. This can result in a catas- 
trophic growth of the mass-transfer and disruption of 
the disk in just a few dynamical timescales. 

Abramowicz et al. ][60l studied the properties of 
mass transfer using many simplifying assumptions: a 
pseudo-Newtonian potential for BH gravity [67J, con- 
stant specific angular momentum of the disk, and ap- 
proximate treatment of the disk self -gravity. They found 
that the runaway instability occurs for a large range of 
parameters, such as the disk radius and the disk-to-BH 
mass ratio Md/Mbh- Subsequent and somewhat more 
refined studies found that the rotation of the BH has a 
stabilizing effect [68 69 1, and a non-constant radial pro- 
file angular momentum was found to strongly disfavor 
the instability [69-72 J. On the other hand, studies us- 
ing a Newtonian pseudopotential for the BH l73l l74l 
and relativistic calculations with fixed spacetime back- 
ground Il65l 1751 found indications of the self -gravity of 
the disk to favor the instability. However, Montero et 
al. (76 J recently performed the first self -consistent and 
fully general relativistic simulations of thick accretion 
disks around BHs in axisymmetry for a few dynamical 
timescales. They found no signatures of a runaway in- 
stability during the simulated time, perhaps implying 
that the self-gravity of the disk does not play a critical 
role in favoring the instability, at least during the first 
few dynamical timescales. 

The problem of the existence and development of 
non-axisymmetric instabilities has a long history. For 
thin Keplerian self-gravitating disks it was found that 
the Toomre parameter [77 - 79 1 can be used to determine 
stability against both local clumping or fragmentation, 



and formation of global non-axisymmetric modes. For 
thick pressure-supported disks, Papaloizou and Pringle 
discovered [61] the existence of a global hydrodynam- 
ical instability that develops on a dynamical timescale 
in disks with negligible self-gravity and constant spe- 
cific angular momentum. A follow-up publication 163 
found this instability also in the disks with power 
law distribution of specific angular momentum £(r) = 
h{rlr c f-1 for all q > V3. Kojima EHSU] found the 
Papaloizou-Pringle (PP) instability in equilibrium tori 
on a fixed Schwarzschild background llBTI using a lin- 
earized perturbative approach. 

Subsequent works clarified the nature of the PP insta- 
bility [82-85J, established how it redistributes specific 
angular momentum [86 1 and discovered that accretion 
has a stabilizing effect on the disk 1 87 , 88 1 . In particu- 
lar, Narayan et al. [82J showed that the PP modes are 
formed by two boundary wave-like perturbations with 
energy and angular momentum of opposite signs that 
are coupled across a forbidden region near the mode 
corotation radius. For wide disks around BHs, the ac- 
cretion suppresses the development of the inner bound- 
ary wave and therefore has a stabilizing effect on PP 
modes [87, 88]. For slender disks, the development of 
PP modes is mostly unaffected by accretion Il88l . The PP 
instability itself amplifies accretion by exerting torques 
on the disk and redistributing specific angular momen- 
tum ISO. 

When the self-gravity of the disk has been taken into 
account, it was found [89-91] that two new types of 
non-axisymmetric instabilities appear, while the PP in- 
stability disappears for most of the models except slen- 
der ones with weak self-gravity. The first type of un- 
stable modes (J-modes) appears in strong self-gravity 
regime, and it is an analog of the classical Jeans insta- 
bility. The second type of unstable modes was found 
in the strong and medium self -gravity regimes [89 1. The 
modes of this type are referred to as intermediate modes 
(I-modes) and represent elliptic deformations of the disk 
(or triangular, square, etc. deformations for higher az- 
imuthal numbers - see e.g. [91 j). 

Yet another type of instability, the so-called "eccen- 
tric instability", was discovered in [92J for thin nearly 
Keplerian self-gravitating disks when the central mass 
was allowed to move. An elaborate mechanism called 
"SLING amplification" [93] was proposed to describe 
this instability. Subsequent investigation 1 94J of this in- 
stability in thin disks suggested a different mechanism 
and predicted that the system is dynamically unstable 
only when the mass of the disk exceeds the mass of the 
central object. 

Finally, Woodward, Tohline, and Hachisu [64 1 pre- 
sented an extensive parameter study of thick self- 
gravitating disks in Newtonian gravity to determine 
the types, growth rates and pattern speeds of non- 
axisymmetric modes. The central mass in their simu- 
lations was allowed to move, and they used 3D time 
evolutions of the disk models with wide range of pa- 



3 



rameters, including disk-to-central object mass ratios 
M D /M C = 0, 1/5, 1 and oo. 

Several recent publications address accretion disks 
and instabilities in these disks in context of GRB cen- 
tral engines. In |20|, Rezzolla et al. studied the proper- 
ties of accretion disks resulting from binary NS merg- 
ers. They obtained thin accretion disks with masses 
~ 0.01 — 0.2 Mq and no evidence of growing non- 
axisymmetric modes or runaway instability. In [95 ], Tay- 
lor, Miller, and Podsiadlowski used 3D SPH simulations 
with detailed microphysics and neutrino transport to 
follow a collapse of an iron core up to the formation of a 
thin massive accretion disk and development of global 
non-axisymmetric modes. They found that torques cre- 
ated by the non-axisymmetric modes provide the main 
mechanism for angular momentum transport, leading 
to high accretion rates of <~ 0.1 — 1 Mq/s, which may 
create a favorable conditions for powering GRBs. 

The stability of disks to runaway and non- 
axisymmetric instabilities is a three-dimensional 
problem that has to be addressed in the framework of 
full GR. One issue of importance is to understand if 
the runaway instability is affected by non-axisymmetric 
instabilities, and vice versa. Despite significant theo- 
retical and computational effort, previous studies of 
the runaway instability do not give a definite answer 
to this question. To our knowledge, only the work by 
Rezzolla et al. [20 1 explored the stability of the disks in 
3D and full GR (for the disks that form in their binary 
NS merger simulations). However, it is important to 
create a comprehensive overall picture of the stability 
of accretion disks for a richer variety of disk models, 
which would require investigating a wider range of 
parameters. In our work, we study in detail the stability 
of slender and moderately slender disks with constant 
distribution of specific angular momentum, which are 
expected to be more unstable both to runaway Il70l ITTl 
and to non-axisymmetric instabilities [87 , 88 ] compared 
to models of Rezzolla et al. [20 1. Our study is based on 
three-dimensional hydrodynamics simulations in full 
GR. 

We model our disks using the ideal fluid approxi- 
mation (i.e. without viscosity) with a T-law equation of 
state (EOS) and T = 4/3. We do not include additional 
physics such as magnetic fields or neutrino and radia- 
tion transport due to the complexity and computational 
cost of the resulting problem. Nevertheless, the adopted 
approach will allow us to identify GR effects which can 
operate in more complex setting that include more re- 
alistic microphysics, neutrino / radiation transport and 
magnetic fields. Also, we limit ourselves to the case of 
non-rotating BHs, while the case of rotating BHs will be 
studied in future publications. 

Another aim of this paper is to estimate the detectabil- 
ity of the gravitational waves (GWs) by the accretion 
disk dynamics. The radial and /or non-axisymmetric 
oscillations of accretion disks can be a source of strong 
GWs. If non-axisymmetric deformations persist for long 



enough time, then the emitted GWs can be detectable 
by current and future GW detectors [96J, provided the 
source is not too far away. Further work on the stabil- 
ity properties of accretion disk could shed more light on 
the number of cycles a non-axisymmetric deformation 
can persist, and thus on the prospects of detecting GWs 
from such systems. 

This paper is organized as follows: Section |n] de- 
scribes the formulations and numerical methods used in 



this paper, including multiblock grids (subsection |II A I 
and the formulations used to evolve the general rel- 
ativistic (subsection II C) and hydrodynamic (subsec- 
tion |II B| equations. Section III describes the grid setup, 
initial data construction procedure and the analysis 
techniques for the non-axisymmetric instabilities. Sec- 
tions IV and [V] present the results of the time evolu- 
tion of self-gravitating disks and the analysis of non- 
axisymmetric instabilities. 

Throughout the paper we use CGS and geometrized 
units G = c = 1. 



II. NUMERICAL METHODS 

The numerical time-evolution scheme used in our 
study can be split into two main parts: the spacetime 
evolution, and the fluid dynamics equations. These two 
parts are evolved in a coupled manner. The numer- 
ical code that we use has been developed within the 
CACTUS computational infrastructure [97 98 1, and uses 
the Carpet mesh refinement and multiblock driver J99j 
I100I1 . A separate module based on Carpet provides a 
range of multiblock systems to represent a variety of 
computational domains for 3D evolution codes |101|. 
The module for evolving the GR hydrodynamics equa- 
tions uses the THOR multiblock code 111021 , which has 
been coupled to the multiblock-based module Quilt 
for evolving the spacetime |103|. The latter implements 
the Generalized Harmonic formulation of the Einstein 
equations in first-order form [104 ] . Below we give a brief 
description of the methods implemented in each mod- 
ule. 



A. Multiblock approach 

For mesh generation and parallelization purposes, 
we employ the multiblock infrastructure developed 
by 11011 . The multiblock approach is widely used 
in astrophysical and numerical relativity simulations 
(see l36lll05rlT23l and references therein). 

Figure [T] shows meridional cuts (cuts in the xz-plane) 
of the block systems that we use. The first two block 
systems on Fig. [T] are the seven-block (also known as 
a "cubed sphere^l24]) and the thirteen-block systems. 
They both have a spherical outer boundary and no coor- 
dinate singularities. The thirteen-block system addition- 
ally has spherical grid coordinate surfaces in the outer 
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cal coordinates. Also, multiblock systems such as dis- 
played on Fig. [Hallow at low cost to extend the computa- 
tional domain far outwards in order to e.g. causally dis- 
connect the system dynamics from the outer boundary 
and to allow more accurate extraction of gravitational 
waves I103U127I . 



B. Hydrodynamics evolution 



(a) 



(b) 





(c) 



(d) 



Figure 1: Cross-sections of the patch systems used in this work 
in the xz-plane: (a) seven-block system, (b) thirteen-block sys- 
tem, (c) six-block system with straight coordinate lines, (d) six- 
block system with distorted radial coordinate lines, described 
in Section lTlI Cl below. 



layer of blocks, which is very convenient for computing 
spherical harmonics and extracting gravitational waves. 
The thirteen-block system is well adapted for simulat- 
ing processes that involve a small source of waves, sur- 
rounded by an extended spherical wave propagation re- 
gion. The central cubical block can have a high resolu- 
tion and can be used to model accurately the dynamics 
of the source, the intermediate layer can represent a near 
zone, while the outer layer of blocks model the radia- 
tion zone, which can carry an outgoing radiation using, 
e.g., constant angular and radial resolutions. These two 
systems are used in the Appendix for test evolutions of 
rotating and non-rotating polytropic stars. 

The remaining block systems on Fig. [TJ; and [TJi are 
six-block systems, and these systems are more suitable 
for modeling configurations with a central black hole. 
The block system displayed on Fig. [TJi is similar to the 
block systems that we use for simulating accretion disks 
around black holes (see Section [ill Q . These two block 
systems contain empty spherical region, which is used 
for excising interior of the BH. 

The multiblock approach has a number of additional 
advantages, including a smooth excision boundary [115 
11251 11261 , a smooth spherical outer boundary suitable 
for radiating boundary conditions, a spherical wave 
propagation zone, and absence of coordinate singular- 
ities such as those associated with spherical or cylindri- 
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The evolution equations of a relativistic fluid are 
derived from the covariant equations of conservation 
of rest mass \7 a (pu a ) — and energy-momentum 
V{,T afo = 0, where T ab has the following form for an 
ideal fluid EH BIO : 

T ab = (p + u + P)u a u b + Pg ab , 

where P is the fluid pressure, p is the rest mass density, 
u is the internal energy density 3 , u" the fluid's four- 
velocity and g ab is the spacetime metric in contravari- 
ant form. The quantities {p, u 1 , u} form a set of primitive 
variables that uniquely determine the state of a single- 
component relativistic fluid at every point in space. 
Our evolution equations read: 

0, 

where f is the time coordinate, g = det(g fl ;,) is the 
determinant of the spacetime metric and T d ac are the 
Christoffel symbols associated with this metric. Af- 
ter introducing a set of conserved variables {D = 
\f-gpu 1 , Q a = \/—gTa}, the equations can be cast into a 
flux-conservative form (as in 11131111321 ; see details spe- 
cific to our particular scheme in [102]): 

3 f D + 3,-D'' = 0, 

d t Q a +d i Q i a = S a , 

where D 1 and Q l a are the fluxes of the conserved vari- 
ables D and Q a , while S a are the source terms for Q a . 

These equations are solved on each block using a fi- 
nite volume cell-centered scheme in the local coordinate 
basis of that block. The reconstruction of the primitive 
variables on the cell interfaces is performed using the 
piecewise-parabolic monotonous (PPM) method 111331 
113411 , while the fluxes through the interfaces are cal- 
culated using a Harten-Lax-van Leer (HLL) Riemann 
solver [135J. In order to compute fluxes, source terms, 
and the stress-energy tensor T ab , the conservative vari- 
ables need to be converted into primitive ones at every 



3 It is also common to define a specific internal energy £ 

[1301 
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timestep. This is done using a Newton-Rhapson iter- 
ative 2Dyy solver [136] with the non-isentropic T-law 
EOS. If at a particular cell on the grid the procedure 
of primitive variables recovery fails to produce physi- 
cally meaningful values (which can happen for number 
of reasons; see, e.g., 111361 1, we use a lDp solver with an 
isentropic polytropic EOS. The vacuum region outside 
the disk is approximated by a low-density artificial at- 
mosphere, whose density is chosen to be 10 -6 — 10~ 7 of 
the initial maximum density in the system. 

The boundary conditions for the hydrodynamics vari- 
ables are imposed on the interblock boundaries using 
interpolation from neighboring blocks. The overlap re- 
gions where interpolation is performed are created by 
adding extra layers of grid points on each block face. At 
the outer boundary and the inner excision boundary, we 
impose outflow boundary conditions. Further details on 
this numerical method of combining multiple blocks can 
be found in 11021. 



lated are absent in the spacetime grid blocks. The ad- 
vantages of using multiple blocks were outlined in the 
Section II A above. The grids on any two neighboring 
blocks in the system are designed to share a 2D inter- 
face grid, where so-called penalty boundary conditions 
are imposed [143 J. For the spatial numerical differen- 
tiation, we employ high-order convergent finite differ- 
encing operators that satisfy summation by parts (SBP) 
property [144]. The SBP property together with penalty 
boundary conditions guarantees strict linear numerical 
stability II143II . All simulations of self -gravitating disks 
with dynamical treatment of relativistic gravity in the 
paper were performed with the finite differencing oper- 
ators of 8th order convergence in the bulk of the grid, 
and 4th order convergence at the boundaries. Time in- 
tegration is performed with a 3rd order accurate Runge- 
Kutta method, that satisfies a total variation diminishing 
(TVD) property 045|. 



C. Spacetime evolution 

To evolve the spacetime metric g a fj we use a gen- 
eralized harmonic (GH) formulation of the Einstein 
equations in the first-order representation developed 
by Lindblom et al. 111041 . In this formulation, the co- 
ordinate conditions are specified using a set of four 
gauge source functions H a , which need to be pre- 
scribed a priori. In our simulations, we have chosen 
the so-called "stationary gauge" [104, 137], in which the 
gauge source functions stay frozen at their initial val- 
ues, H a (t, x l ) = H a (to, x l ). Such a choice of the gauge is 
convenient for quasi-stationary spacetimes, such as per- 
turbed BHs 1 103 1 or accretion disks around BHs. 

The first-order representation that we use is linearly 
degenerate, symmetric hyperbolic [138-140], and con- 
straint damping. Linear degeneracy guarantees that the 
system will not develop gauge shocks during the evolu- 
tion 11411 . Symmetric hyperbolicity with boundary con- 
ditions imposed on incoming characteristic fields is the 
necessary condition for well-posedness of initial bound- 
ary value problems (see, e.g., Chapter 6 of Gustafsson et 
al. USD- 
Being constraint damping, which is one of the most 
important advantages of the chosen formulation, means 
that the constraints on the variables are included into 
the evolution equations in such a way that they are ex- 
ponentially damped during evolution. Such a property 
not only reduces the error in the solution, but also elimi- 
nates a variety of numerical instabilities associated with 
unbounded growth of the constraints. 

The spacetime metric evolution equations are dis- 
cretized using finite differences on multiple grid blocks 
that cover the computational domain 111011 11021 . This 
means that the spacetime grid blocks are only sub- 
sets of the hydrodynamics grid blocks; the overlap re- 
gions where the hydrodynamics variables are interpo- 



D. Code tests 

Numerical methods for solving the general relativis- 
tic hydrodynamics equations and the spacetime equa- 
tions are inherently complex, and codes need to be thor- 
oughly tested before they can be successfully applied to 
physical problems (see e.g. 11461 114710 . The hydrody- 
namics code THOR was tested for fixed spacetimes by 
Zink et al. in |102|, where it was demonstrated that the 
code can handle situations encountered in many astro- 
physical scenarios, including relativistic shocks, rotating 
polytropes and equilibrium tori around BHs 1811 . The 
spacetime evolution code Quilt was tested in [ 103"l ll48l . 

We present results of our tests of the coupled hy- 
drodynamics and spacetime evolution codes in Appen- 
dices ^][3 In particular, in Appendix [2] we report sta- 
ble convergent evolutions of a rapidly rotating poly- 
tropic star, and in Appendix [3| we demonstrate that our 
code faithfully reproduces fundamental frequencies of a 
non-rotating polytropic star. Additionally, recently Zink 
et al. [149] have used the coupled THOR+QuiLT code 
to measure the frequencies of /-modes and to study 
their neutral points in the context of the Chandrasekhar- 
Friedman-Schutz instability. 



E. Constraint damping 

The constraint damping scheme for the spacetime 
evolution mentioned above depends on two freely- 
specifiable parameters, k and 72. These parameters can 
be freely chosen as a function of space and time, i.e. the 
system of equations does not directly depend on their 
spatial or temporal derivatives. Here we describe how 
these parameters are specified in our simulations. 

In the generalized harmonic formulation, the set u a = 
{gab'^iab'^ab} or dynamical fields consists of the met- 
ric g a f, and linear combinations of its derivatives <3> !a {, = 
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digab> ^ab — ~ n<: ^cgab> where n c is a normal to a t — 
const, hypersurface. The evolution equations for the 
fields can be written in the following form: 

d tgab = G ab (x a ,u a ) r (1) 
d t <S> iab = T iab (x a , u a ) - j 2 C iab , (2) 
d t U ab = P ah (x a , u a ) - KC ab - JzPQab, (3) 

where C ab , G{ ab are the constraints, S' is the shift vec- 
tor, and G ab , T{ ab , P ab are right-hand sides of the for- 
mulation without the constraint damping terms. In the 
continuum limit, the constraints are zero, and the sys- 
tem iJTJ reduces to the original Einstein equations. At 
the discrete level, the constraints can be non-zero, and 
the constraint damping terms provide a non-vanishing 
contribution to the right-hand sides of the system fiT. 
Since constraint damping terms should act only as small 
corrections to the evolution system, their contribution 
should not exceed those of functions G ab , 7-^, P ab - Oth- 
erwise, the evolution of the system will be dominated 
by the numerical constraint violations. 

We notice in our simulations that the functions T{ ab , 
P ab on the right-hand sides of ([lj fall off as l/r a with 
a ~ 4 beyond r > r^isk' where is the approximate 
outer radius of the disk. At the same time, the constraint 
violations C ab , Ci ab fall off as l/r" with B ~ 1. This 
means that, if we use constant values for 72 and k, the 
constraint damping terms will dominate the dynamics 
of the system for sufficiently large r. In order to avoid 
this situation, the functions k and 72 must fall off with 
the radius as ~ 1/r 3 . 

We find that the following radial profiles of K and 72 
for our simulations lead to satisfactory results: 



K(r) = 7C* 

72 0) = 7* 



7r v 1 + n 



+ arctan r* 



+ arctan r* 



(4) 
(5) 



where r* = (r — ro)/o~, and k*, 7*, yq, <t are (posi- 
tive) constants. This radial profile approaches a constant 
value of ~ k* (or ~ 7*) for r < tq, and falls off as ~ 1/r 3 
for r ^ tq. The parameter u determines the extent of 
the smooth transition region between these two regimes. 
We use = 7* = 4, r$ = 12, and a = 8 in our simu- 
lations. Such profiles of the constraint damping coeffi- 
cients allows imposing strong constraint damping near 
the BH and the disk without introducing spurious dy- 
namics far away from the origin. 



III. INITIAL SETUP 

We set up initial equilibrium disk configurations by 
solving Einstein constraints using a version of the RNS 
solver 1 150 1 adapted to the problem of equilibrium 
tori. The method of solution is similar to the one used 




Figure 2: Conformal Penrose diagram of an axisymmetric 
spacetime, consisting of a non-rotating BH, distorted by a mas- 
sive stationary disk around it. Each point on the diagram cor- 
responds to a spheroid, located at a given geodesic distance 
from the BH horizon. Blue lines represent a quasi-isotropic 
foliation of the spacetime, while red lines show a horizon- 
penetrating foliation. 



in [151]. The spacetime is assumed to be stationary, ax- 
isymmetric, asymptotically flat and symmetric with re- 
spect to reflections in equatorial plane. The metric is a 
general axisymmetric metric in quasi-isotropic coordi- 
nates: 



ds 2 



■■ -X 2 dt 2 + e 2K {drl + r 2 je 2 ) (6) 
B 2 1 X 2 r\ sin 2 6 (df-cvdt) 2 . (7) 



where (t,r*,8,cp) are the coordinates, and A, a, B and 
co are metric potentials which depend only on and 
9. The RNS code implements the KEH(SF) method 11501 
115211 , in which the Einstein equations for the metric po- 
tentials A, B and to are transformed into integral equa- 
tions 1 151 1 using Green's functions for the elliptical dif- 
ferential operators. The remaining metric potential a 
can then be found by integrating an ordinary differential 
equation, once the rest of the potentials are known. The 
KEH(SF) method uses compactified radial coordinate s 
which maps the region [0, 00) into a segment [0,1]: 



s = 

r + r + 

where r+ is the outer radius of the disk. The boundary 
conditions are imposed at symmetry interfaces and at 
the event horizon. The latter can always be transformed 
to a sphere of compactified radius S/,, while preserving 
the form of the metric given above 11511 . At the horizon, 
we impose the boundary conditions for a non-rotating 
BH: A = B = co = 0, and we set B/A = e a at the sym- 
metry axis. The corresponding integral equations are 
then solved using Newton-Raphson iterations [150 1 in 
the upper quadrant 6 £ [0, n/2], s £ [0, 1] of the merid- 
ional plane. 

The resulting quasi-isotropic metric is degenerate at 
the event horizon, which is very problematic for the evo- 
lution with excision of the BH interior using the gener- 
alized harmonic formulation, since this method requires 
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coordinates without pathologies at the event horizon. 
This situation is best illustrated by a conformal picture 
of the complete spacetime, shown in Fig. [2] All quasi- 
isotropic slices meet the horizon of the BFiat its throat, 
which makes the metric degenerate at the horizon. Re- 
gions II and IV are not covered by the quasi-isotropic fo- 
liation. It is necessary for our time evolution methods to 
have a time-independent foliation which penetrates the 
horizon and continues in region II rather than region III. 
There are several options to address this issue: 

(1) solve the complete system of equations in 
horizon-penetrating coordinates rather than 
quasi-isotropic ones; 

(2) use puncture initial data, as developed in 1 153 J, 
and choose such a gauge for the evolution that af- 
ter some time the spatial slices move from region 
III to region II, similarly to what happens with 
punctures in the BSSN formulation and 1 + log 
slicing 111541 ; 

(3) apply a spacetime coordinate transformation 
from quasi-isotropic to horizon-penetrating coor- 
dinates. Since there is still no solution provided 
in region II, it needs to be extrapolated into that 
region. 

We use option (3), since it is easier to implement in 
the context of the generalized harmonic system, and fits 
more naturally in the gauge choice employed by RNS. 
Details on the spacetime transformation that we apply 
to the initial data can be found in Appendix [T] 

A. Blending numerical and analytical metrics 

We could not extrapolate the initial data produced by 
the elliptic solver to the region II inside the horizon, be- 
cause the data does not have enough smoothness near 
the horizon. The problem with initial data at the horizon 
seems to be similar to the problem with Gibbs phenom- 
ena, when solutions exhibit an oscillatory behavior and 
a lower order of convergence near a stellar surface. To 
handle this problem, we use an approximation, in which 
we replace the numerical metric in the region where it is 
not accurate with an analytic Kerr-Schild solution of the 
same BH mass. We blend the numerical metric g'^ m and 
the Kerr-Schild metric g^ b s using the following prescrip- 
tion: 

g ab = (l-w(r))g n a r + ™(r)gi S 
where the weight function w(r) is defined as: 

fl, if r<b lr 

w{r) = I cos 2 2(b~-b!) ' if r e H. 
|o, itr>b 2 , 



and the segment r 6 [^1,^2] determines a finite-size 
blending zone between the two metrics. The weight 
function w{r) is non-constant only in a narrow spheri- 
cal layer outside the horizon. 

Blending two metrics introduces constraint violations 
at the continuum level, but they subside rapidly in time 
due to the constraint damping property of our evolution 
scheme. The constraint damping scheme, however, does 
not necessarily satisfy the conservation of mass, so after 
the constraint violations are suppressed, the system ar- 
rives at a different state, which can be characterized as a 
close equilibrium configuration with some axisymmet- 
ric gravitational perturbation. 

The location and size of the blending layer can be ad- 
justed to minimize the initial unphysical oscillation in a 
BH mass. For the evolutions presented below, we used 
blending in the range r £ [1.05 r g , 1.15 r g ], where r g 
is the radius of the BH event horizon. Such a choice 
results in the initial oscillation of the BH mass within 
« 12%, which then settles down to a stationary value 
of 97.5 ± 0.5% during the first orbital period (see Fig.|8|3 
and related discussion in Section [TV). 



B. Initial disk models 






Figure 3: Contours of the disk surfaces for models A, B and C, 
and the location of BH horizon in the meridional plane. 



We have constructed three initial disk models with 
disk-to-BH mass ratios M d /M B h = 0.235, 0.174 and 
0.108, labeled A, B and C. Model C is slender, which 
means that its width is much smaller than the radius 
of the torus, i.e. r c 3> (r+ — r_). Models A and B are 
moderately slender, i.e. for these models r c w (r+ — r_). 
Figure [3] shows the contours of the disk surfaces for our 
models and the location of BH horizon in the meridional 
plane. All models are constructed using a polytropic 
EOS with polytropic index T = 4/3, constant specific 
entropy, and a constant specific angular momentum dis- 
tribution. 
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Model 


A 


B 


c 


Specific angular momentum £ [Algjf ] 


4.50 


4.32 


4.87 


T-*nl \rtrnni r rnnQfant flO^^rm^cr 2] 

1 U1V UUplL I^UILDICIILL IV 1U * 111 J 


1.28 


1.04 


0.519 


MaYiiniim Hpirsitv n M D^^cr /rtn 1 

IVldAUl 1 Llll I LlCllDlly kJ£ ±\J i \ 111 


1.23 


1.17 


0.755 


T~)i q V-tn-Rl-T maiQ ratio. A/Tn / AAmi 


0.235 


0.174 


0.108 


Kinetic to potential energy T / 1 W| 


0.479 


0.497 


0.499 


Central radius Tq I Tg 


6.51 


5.58 


8.23 


Ratio of inner to central radius r_ / r c 


0.655 


0.655 


0.756 


Ratio of inner to outer radius r_ / r+ 


0.385 


0.407 


0.600 


Orbital frequency at r c , O c [s _1 ] 


1713 


1912 


1121 


Self-gravity parameter T = inGp/Ci^: 


3.52 


2.68 


5.04 



Table I: Physical parameters of the self-gravitating initial disk 
models, used in our simulations. 



Table [H lists physical and geometrical parameters 
of the disk, including the ratio of disk-to-BH mass 
Md/Mbh/ the ratio of kinetic energy T to potential en- 
ergy W and the self-gravity parameter t. The latter can 
be defined as p c /p S piu where p s ^ = Q^/47rG is the 
density of a uniform sphere with radius r c that creates 
equivalent gravity at that radius. Notice that T/ \ W | cor- 
relates to the "slenderness" of the disk r_ /r+ . The value 
of specific angular momentum I is given in units of the 
BH mass. These parameters allow us to make qualita- 
tive and quantitative comparisons between our mod- 
els and the models studied in previous works l64l l89l - 
[9TirT55l. 



C. Adapted curvilinear grid 

In order to accurately resolve both the disk and the BH 
while minimizing the computational cost, we have de- 
signed a series of curvilinear multiblock grids adapted 
to each of the disk models. To obtain such a grid, we 
start from a six-block system (displayed on Fig Jit) that 
was previously used in QUILT for computationally effi- 
cient and numerically accurate simulations of perturbed 
BHs 1 103 1 . We apply a radial stretching and an angu- 
lar distortion to the six-block system so as to create a 
uniformly high resolution near the BH, nearly cylindri- 
cal grid near the disk, and approach a regular six-block 
spherical grid in the wave zone. These mappings are 
described in detail below. Figure |TJi shows an illustra- 
tion of the grid distortion, while Fig. [4^ shows the actual 
curvilinear grid used in some of our simulations. 

A regular six-block system consists of two polar 
blocks (near the z-axis) and four equatorial blocks (near 
the xy-plane). We can assign quasi-spherical coordi- 
nates {r,6,cp} and {?", 0i, 02} to the equatorial and polar 
blocks, respectively. They can be related to the Cartesian 
coordinates {x,y,z} by the following transformation: 

(a) for an equatorial block in the neighborhood of the 



positive x axis: 

x = rl sf\ + tan (p 2 + tan 2 , 
y = x tan cp, 
z = xtan0, 

(b) while for a polar block in the neighborhood of the 
positive z axis: 

x = z tan 0i, 
y = z tan 02, 

z = rf sj\ + tan 0i 2 + tan 2 2 . 

(c) The remaining blocks are obtained by applying 
symmetry transformations. 

We set the coordinate ranges for the polar blocks to be 
r e [Rmin,Rmax],0i,Q2 & [—0*, 0*], where the value of 0* 
controls an opening angle of the polar blocks. For equa- 
torial blocks, r e [R min , R max ], 8 G [-§ + 0*, f - 0*] 
and cp E [— -j , j ]. The size of the numerical grid for 
each block is fixed by three numbers: N r , Nq, and Ng. 
The equatorial blocks have N r x Nq, x Nq cells, and the 
polar ones have N r x Ng x Ng cells. 

In order to obtain a variable radial resolution, we 
apply a smooth one-dimensional radial stretching S : 
r f that yields the desired resolution profile A r (r). 
This profile is chosen based on several stringent require- 
ments imposed by an accuracy and available computa- 
tional resources. First, there need to be at least 8 grid 
points between the excision radius and the BH horizon 
in order to prevent constraint violations from leaving 
the interior of the BH (this question is considered in 
some detail in [156]), as well as to allow for some (re- 
stricted) BH movement. Second, our convergence tests 
show that near the disk, the grid needs to allow for at 
least 40 points across the disk in order to achieve a global 
convergent regime in hydrodynamical evolutions. For 
this purpose, the radial resolution profile is adapted to 
have sufficiently high resolution near the disk as well. 
Third, in the wave zone, there is no need to maintain 
very high radial resolution. However, this resolution 
needs to be uniform (rather than, for example, exponen- 
tially decreasing) to be able to carry the radiation accu- 
rately without dissipation. All these requirements result 
in the radial resolution profile A r (r) shown on Fig.|Ij). 

The resolution in the 8 direction Ag near the disk also 
needs to have at least 40 points across the disk for con- 
vergence. However, simple increase of Ng in the six- 
block system leads to a very small minimal grid step 
A m! „ at the excision radius R m i n that is too restrictive for 
the time step due to CFL condition. To increase the grid 
resolution in the 0-direction near the disk only, we apply 
a radial-dependent distortion to the angular coordinates 
0,01 and 2 : 

0(0» =0(l-/3(r)f^), 

0,(0--, r) =0-(l + j6 (r)^), i = l,2 
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where the function /3(r) is the amplitude of the dis- 
tortion, chosen to have a Gaussian profile /5(r) = 



^* exp( 



tq) 1 1 cr ), in which the parameters /3*, tq 



and a are chosen to satisfy the above requirements. This 
distortion bends diverging radial coordinate lines to- 
wards equatorial plane around the radius Tq, making the 
grid resemble a cylindrical shape near that radius (as il- 
lustrated by Fig. Efe). Figure Wp shows the dependence 
of A# on r along the x-axis (dotted black line) and along 
the z-axis (dash-dotted red line). The radial-dependent 
^-distortion increases Ag on the x-axis at the expense of 
Ag on the z-axis near the disk. Away from the disk, the 
distortion vanishes and Ag approaches the linear depen- 
dence Ag(r) a r. 

Figure|4 shows an example of the resulting curvilinear 
grid. For this example, 0* = 30°, Ng = 49 and N f = 25. 
These parameters make the grid spacing approximately 
uniform in angular direction near the BH and at large r. 
The inner (excision) radius is R ml >, = 1.3 and the outer 
one is R m ax = 500, while the minimum grid spacing is 
A min = Ag(r m! „) 0.06. The apparent horizon has ra- 
dius of f» 1.7, so that as many as 8 grid points can be 
placed inside the horizon without decreasing the mini- 
mum grid spacing. Distortion parameters of the Gaus- 
sian have values Tq — 12, a = 8, and /3* = 1.3, which 
is sufficient to concentrate about 40 grid points across 
the disk in the vertical direction. The radial resolution 
is adapted to be (A t )bh « 0.06 at the excision sphere, 
(Ar)disk ~ 0-25 around the disk, and (A,-)wz ~ 5.2 in the 
wave zone. 

All of the time evolutions described in the next two 
Sections [IV] and [V] use adapted six-block curvilinear 
grids. Table [TT] summarizes dimensions and resolutions 
of the grids for each of the simulations used in these 
two sections. The first column lists simulation names, 
which consist of two or three symbols. Models K1-K6 
are considered in Section IV Al In the rest of the model 
names, the first letter denotes the initial disk model (A, 
B or C), the second letter signifies whether the simula- 
tion is evolved in full GR (F) or Cowling (C) approxi- 
mation (see corresponding Sections VC and |V B| . The 
next number (if present) is the azimuthal number m of 
an added non-axisymmetric perturbation (as explained 
in Section IV). Simulations AFc, AF and AFf represent 
disk model~A evolved in full GR using coarse, medium 
and fine resolution grids. These models are used in Sec- 
tor convergence studies. 



tion 



IV 



Finally, note that the 
adapted curvilinear grid example above (displayed in 
Fig. 4} corresponds to the the grid used in simulation AF, 
and the grids of simulations AFf and AFc are obtained 
by increasing and decreasing the resolutions in the AF 
grid by a factor of 3/2. 



linear perturbative studies of accretion disks (e.g. in 
l80l ). Namely, we analyze first few terms in Fourier ex 
pansion in angle cp of the disk density p(r,cp) on a se 
quence of concentric circles in the equatorial plane: 



P(t,r,<p)=p{t,r) 1 + 



m—\ 



-i(w m t—m<p) 



where p is a <p-averaged density at a given radius. The 
quantity D m represents the (complex) amplitude of an 
azimuthal mode m, the real part of the quantity co m de- 
termines a mode pattern speed, while its imaginary part 
determines a mode growth rate. Following |64, 153, we 
quantify the growth rate and the pattern speed of a non- 
axisymmetric mode by two dimensionless parameters 
j/i and 1/2/ defined as 



y t (m) 



Re(w m ) 



a 



orb 



m, 



3/2 (m) 



n 



orb 



We calculate a value of the parameter )jx{m) from a slope 
of log | D m | versus t line at an arbitrary radius, while 
y\{m) is obtained from a slope of the mode phase an- 
gle (pm = (pm{i)- Notice that because the modes that 
we consider are global, their growth rates and pattern 
speeds do not depend on a radius. 

We use the parameter \j\{jn) to obtain the value of a 
corotation radius r cr for a given mode: using the mode 
pattern speed, 



Qp — d rb 



yiO) 



m 



and the radial profile O — Ci(r) of a fluid angular ve- 
locity in the disk, we can calculate r cr from the equation 
Cl(r = r cr ) = n p . 

In those simulations where the disk is oscillating ra- 
dially, the values of the mode amplitude D m at a given 
radius oscillate due to disk oscillations, which makes it 
more difficult to extract the mode growth rates from D m . 
In such cases, we have found that more accurate growth 
rates can be obtained if we use normalized root mean 
squared (RMS) mode amplitudes G m , defined as: 

G m = (D m ) 2 I (Do) 2 

denote an RMS value 



where the angle brackets (. . . ) 2 
over radii from r_ to r+: 



(D 



1 



ml 2 



D m \ 2 dr 



1/2 



D. Data analysis 

To identify and characterize non-axisymmetric insta- 
bilities, we will adopt an approach commonly used in 



IV. TIME EVOLUTION 

In this section, we present the results of the fully gen- 
eral relativistic time evolution of the initial data for the 
reference model A, constructed as described above in 
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Figure 4: The left panel (a): meridional cut of an adapted curvilinear grid used in one of our simulations, combined with the 
logarithmic density contours of the disk at t = 0. The intersecting radial coordinate lines belong to the neighboring blocks in the 
overlapping (interpolated) regions. The thick black circle marks a location of an apparent horizon of the BH. The inset in the lower 
left corner of the plot shows the high-resolution grid around the BH, adapted in such a way that the resolution is uniformly high 
in every direction. The right panel (b): radial profile of the resolution in the r (red solid line) and the 9 (dotted lines) directions 
for the grid shown in the top panel. The radial coordinate is plotted in logarithmic scale. The top horizontal axis shows the 
values of the radius in numerical units. The marks on the lower axis are: r ex is the radius of the BH excision sphere; r g is the 
gravitational radius of the BH; r_, r + is the inner and outer radii of the disk; r c is the radius of the density maximum; ryv is the 
"wave extraction" radius; R ou t is the outer radius of the domain. The black dotted line is the resolution in the ^-direction on the 
x axis, and the red dashed-and-dotted line is the resolution in the (9-direction on the z axis. 



Model 


r g 




Rmax 


N x x N z x N r 




{^r)disk 


(Ae)r„* 


(Afl)»ws 


(Ar)wZ 


Kl 


2.0 


4.0 


30 


25 x 49 x 96 


0.07 


0.25 


0.32 


0.72 




K2 


2.0 


4.0 


25 


25 x 49 x 96 


0.07 


0.25 


0.32 


0.72 




K3 


2.0 


4.0 


22 


25 x 49 x 96 


0.09 


0.20 


0.17 


1.28 




K4 


2.0 


6.0 


17 


25 x 49 x 96 


0.08 


0.18 


0.13 


1.68 




K5 


2.0 


7.0 


16 


25 x 49 x 96 


0.03 


0.10 


0.05 


0.50 




K6 


2.0 


8.0 


14 


25 x 49 x 96 


0.02 


0.08 


0.03 


0.53 




AC, AC1, AC2 


1.760 


1.7 


25 


37 x 73 x 144 


0.05 


0.17 


0.21 


0.48 




BC, BC1, BC2 


1.806 


1.7 


25 


37 x 73 x 144 


0.05 


0.17 


0.21 


0.48 




CC, CC1, CC2, CC3 


1.812 


1.7 


25 


37 x 73 x 144 


0.05 


0.17 


0.21 


0.48 




AF, AF1, AF2 


1.760 


1.3 


500 


25 x 49 x 280 


0.06 


0.25 


0.22 


0.95 


5.2 


AFc 


1.760 


1.3 


1000 


17 x 33 x 180 


0.09 


0.37 


0.33 


1.43 


7.8 


AFf 


1.760 


1.3 


1000 


37 x 73 x 420 


0.04 


0.17 


0.15 


0.62 


3.4 


BF, BF1, BF2 


1.806 


1.3 


1000 


25 x 49 x 280 


0.06 


0.25 


0.22 


0.95 


5.2 


CF, CF1, CF2, CF3 


1.812 


1.3 


1000 


25 x 49 x 280 


0.06 


0.22 


0.22 


1.01 


5.2 



Table II: Parameters of numerical grids for all the simulations of accretion disks used in this study. The naming convention of 
the simulations is explained in the main text in Section III C on adapted curvilinear grids. All resolutions and linear dimensions 
are given in computational units. The first column, which lists the values of the BH gravitational radius r ? , allows to convert all 
quantities from computational to CGS units. The remaining columns contain: R m i„, R m ax are radial extents of the computational 
domain; N x is the number of grid points in the horizontal x- or y-direction; N z is the number of grid points in the vertical 
direction; N r is the number of grid points in the radial direction; A ml -„ is the minimal grid step size, in 9 direction at the excision 
sphere; (Ar)disjfc is the resolution in the r-direction at radius r c ; {Ag)r ci x is the resolution in the (9-direction on the x-axis at radius 
r c ; (Ag) rc/Z is the resolution in the 0-direction on the z-axis at radius r c ; (A r ) wz is the radial resolution in the wave zone (not used 
for simulations on fixed background). 
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Section III Overall, the dynamics of model B is qualita- 
tively similar to that of A, while model C exhibits a qual- 
itatively different time evolution. The BH initial mass in 
all of our models is 2.5 -A4q, as in some of the previ- 
ous works [66 , 72, 158], and the disk rotational period 
t c = 3.667 ms is used as a unit of time in all of the plots 
in this section. 

At the beginning of time evolution, the metric blend- 
ing procedure (described earlier in Section III A I intro- 
duces axisymmetric constraint-violating perturbation to 
the spacetime near the BH, causing an unphysical oscil- 
lation in the BH mass, which damps out in about one 
orbital period of the disk. We discard the first orbital 
period when analyzing simulation data and drawing 
physical conclusions about the system dynamics. 

In the meanwhile, the axisymmetric perturbation 
propagates outwards and triggers axisymmetric disk os- 
cillations. Disk oscillations lead to formation of shock 
waves, which transform kinetic energy of the shock to 
thermal energy, resulting in damping of the oscillations. 
We discuss different aspects of the dynamics of the disk 
in more detail below. 

After about three orbital periods, the disk develops an 
m — 1 non-axisymmetric mode, which we have identi- 
fied as Papaloizou-Pringle (PP) type instability [6lJ l62l 
(discussed in more detail below in Section[V}. The same 
type of mode develops in model B, while the more slen- 
der model C develops anm =2 mode of an intermedi- 
ate type (see Section VC[ . As the m = 1 mode grows, 
the center of mass of the disk drifts away from its ini- 
tial position along a spiral-like trajectory, as shown in 
Fig. p] (red squares). As a result of gravitational interac- 
tion between the deformed disk and the BH, the latter 
also starts spiraling away, mirroring the motion of the 
disk center of mass, as shown in Fig. [5] (black line). This 
plot also shows dashed lines that connect the positions 
of the center of mass of the BH and the disk at different 
moments of time. As we can see, all of the these lines in- 
tersect with each other at one point at the initial location 
of the center of mass of the disk-BH system, implying 
that the center of mass of the system does not move, as 
should be the case for BH motion caused by physical in- 
teraction with the disk (but not gauge effects). 

Figure [6] shows the time evolution of the amplitude 
| D\ | of the m = 1 PP non-axisymmetric mode and the 
distance rgn (normalized to r g ) from the BH center to its 
initial position. As we can see, both \D\ \ and rgn have 
the same growth rate. This feature provides another ev- 
idence that the BH motion is a result of the physical in- 
teraction with the m = 1 disk deformation, but not due 
to gauge effects. 

We point out that, as a result of the interaction of 
the m = 1 deformed disk with the BH, the latter ac- 
quires significant orbital angular momentum from the 
disk. Fig. [7] illustrates how angular momentum of the 




x/r n 



Figure 5: Trajectories of the center of mass of the accretion disk 
(red line) and the BH (black line). Dashed lines show four 
consecutive simultaneous locations of the two centers of mass, 
and a small red circle at the origin marks the location of their 
common center of mass. The spiral motion of the BH is caused 
by the development of the non-axisymmetric m = 1 mode in 
the disk. 
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Figure 6: Time evolution of the amplitude |Dj | of the non- 
axisymmetric m = 1 mode in the disk and the length of the BH 
position vector r^n/rg. The growth of these quantities is cor- 
related, i.e. they develop at the same time and with the same 
rate. 
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Figure 7: Illustration of the angular momentum transfer from 
the disk to the BH. The solid (red) line shows the orbital angu- 
lar momentum of the BH /g^ , the dotted (magenta) line shows 
the decrease in the angular momentum of the disk } D — /n, and 
the dashed (blue) line shows the total decrease of the angular 
momentum of disk+BH (/d + Jbh) — Jo- All quantities are di- 
vided by the initial value of angular momentum of the disk Jq. 
The noticeable angular momentum transfer can only be seen 
in the last orbital period, when the m = 1 distortion reaches 
significant amplitude. The gradual decrease in total angular 
momentum is due to the mass loss at the interpolation bound- 
aries. This is a numerical artifact which converges away with 
resolution. 



disk (J D ), BH (/bh) 4 and BH+disk system (J BH + /d) 
changes with time with respect to initial disk angular 
momentum (Jo). The total angular momentum of the 
disk+BH system decreases by ~ 1.5% in ~ 7 orbital pe- 
riods due to numerical errors such as interpolation at 
the block boundaries and evaporation to the artificial at- 
mosphere in the case of medium resolution. As a result 
of angular momentum transfer, Jd / Jq additionally de- 
creases by ~ 1.5%, which is completely compensated by 
the ~ 1.5% increase of /bh/ Jo- 

Unfortunately, the continued outspiraling motion of 
the BH ultimately leads to the intersection of the ap- 
parent horizon with the excision boundary. At this 
point, the inner excision boundary conditions become 
ill-posed, and we have to terminate our simulations. 

Similarly to the findings of l20l l76l , we did not ob- 
serve runaway instability in all three models. We also 
do not expect this instability to occur at a later time. In 
the most likely scenario of the subsequent evolution, the 
development of non-axisymmetric instabilities will re- 
distribute disk angular momentum and lead to a pro- 



file of specific angular momentum that increases out- 
wards Il86l l9ll . Such angular momentum profile was 
shown to strongly disfavor runaway instability II 70 1 . 
Moreover, the damping of radial disk oscillations in 
models A and B reduces and eventually completely ter- 
minates mass transfer from the disk to the BH, prevent- 
ing onset of runaway instability. 
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Figure 8: Time evolution of the constraints and the BH mass 
for simulations with the coarse, medium and fine resolution 
grids, (a) Loo- and Lx-norms of the Hamiltonian constraints 
as a function of time. For the coarse and medium resolution 
cases, the plot shows that the constraints are reduced down to 
the discretization error level during the first orbital period (see 
the main text for more detailed discussion), (b) Time evolution 
of the BH mass (measured from the area of apparent horizon, 
normalized by its initial value). After the initial transitional 
oscillatory phase, the BH mass settles down to 97.5 ± 0.5% of 
its initial value. The oscillation in the first orbital period is 
caused by the constraint violations due to blendi ng pro cedure 
in the construction of our initial data (see Section [lII A| . 



4 We calculate the total angular momentum of the disk using expres- 
sion Jo = J d 3 Xy/—gT' 11591 , while the orbital angular momen- 
tum of the BH is calculated using a simple Newtonian estimate: 
using the BH speed r<p and its distance from the origin r, we get 
Jbh - M BH r 2 (j). 



To determine how much our results depend on nu- 
merical resolution, we have performed simulations of 
this model for three different resolutions with grid cell 
size scaling as 1 : 1.5 : 1.5 2 . We refer to these as 
the coarse, medium and fine resolutions hereafter. The 
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curvilinear geometry of the blocks was chosen to be the 
same for all of the three resolutions. We list the parame- 
ters of the resulting coarse, medium and fine (denoted as 
AFc, AF and AFf, respectively) grid models in Table [TT] 

Figure [8^ shows the L%- and Loo-norms of Hamilto- 
nian constraint violation for the coarse, medium and 
fine resolution simulations. This plot shows that there 
are two distinct regimes in the evolution of the con- 
straints: initial exponential decrease and subsequent 
steady plateau. The latter is due to the constraint damp- 
ing mechanism of our evolution scheme, in which the 
rate of production of discretization errors is balanced 
by the rate of constraint damping. Since the discretiza- 
tion error depends on grid cell size, the value of the 
plateau decreases with increasing resolution. In the case 
with the highest resolution, the "plateau regime" is not 
reached during the time span of the simulation. Fig. |§k 
also shows that there is a small region of a rapid growth 
of the constraints at the very end of the simulations. This 
increase is caused by the approach of the BH apparent 
horizon too close to the inner excision boundary as a re- 
sult of the interaction of the BH with the tn — 1 defor- 
mation of the disk described above. 

Figure shows the time evolution of the BH mass 5 . 
Because the constraint violations due to the metric 
blending are introduced at the continuum limit, the 
BH mass shows an unphysical oscillation that does not 
converge away with resolution, but which completely 
damps out in about one orbital period. The BH mass 
then stabilizes at 97.5 ± 0.5% of its initial value and re- 
mains near this value until the end of the simulation. 
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Figure 9: Time evolution of the maximum rest-mass density 
and the location of a rest-mass density maximum r c = r(p max ). 



We now proceed to the analysis of the axisymmetric 
disk oscillations induced by the metric blending proce- 



dure (described in Section III A) . Fig. [9] shows the time 
evolution of the maximum disk density p max (solid line) 
and the radial position of the disk center r c 6 (dashed 
line), normalized to their values at t = 0. As can be seen 
from this plot, the evolution of both pmax and r c is dom- 
inated by a single mode of oscillation. The disk expands 
while moving away from the BH, and contracts while 
moving towards the BH. This can also be seen in Fig. 10 
which shows meridional cuts of the disk in the xz-plane 
in different phases of the oscillation. 

The radial oscillation of the disk has frequency of ~ 
164 Hz, a value that is slightly smaller than the epicyclic 
frequency k = 201 Hz at r = r c for a Schwarzschild 
BH of the same mass 7 . Rezzolla et al. B162I studied rel- 
ativistic axisymmetric oscillations of accretion disks in 
Cowling approximation using disk models that have the 
same T (= 4/3) and (constant) specific angular momen- 
tum profile as the ones studied here. They found that the 
disks with r c = 3.860 To have the radial oscillation fre- 
quency of ~ 261 Hz (= 0.02017 for their c2 model in nor- 
malized units |162J), which is higher than the ~ 164 Hz 
frequency that we obtain in our case. Part of this differ- 
ence stems from the difference in the disk models: while 
model of 1 162 1 has r c = 3.860 r ? , our disk model A has 
r c = 6.51 r g (cf. TableJIf. Rezzolla et al. Il62l has found 
that the oscillation frequency / depends on the disk ex- 
tent L in a particular way. Namely, / decreases with in- 
creasing L, and in the limit of very thin disks (L — > 0), 
/ approaches the local value of the epicyclic frequency k 
(see Fig. 4 of [162J). If we assume this dependence f(L) 
to hold for any r c , we can obtain the value of radial os- 
cillation frequency for a disk with the same L and r c as 
our model. This frequency turns out to be ~ 168 Hz, 
which is very close to / = 164 Hz that we obtain in our 
simulations. This result might indicate that for the disks 
that are similar to the ones considered here, the disk self- 
gravity does not significantly affect the frequencies of 
the radial oscillations. Note that in the case of NSs, the 
frequency of the fundamental quasi-radial modes was 
found to differ by a factor of ~ 2 between Cowling and 
full GR simulations [146 , 163J, perhaps implying that the 
self-gravity plays a more important role in the case of 
NSs. However, since it is unclear whether such depen- 
dence of / on L holds exactly for other values of r c , this 
result should be taken with caution. We will revisit this 
issue in a future publication. 

The disk oscillations are damped due to formation of 
shock waves that convert the kinetic energy of the os- 
cillations into the thermal one. At the end of the con- 
traction phase of the disk, the high-density inner part 



5 We measure irreducible BH mass Mbh using the area of apparent 
horizon 11601. 



6 We define the radius of the disk center r c to be the radius, where 
the tp-averaged disk density reaches its maximum. 



7 The relativistic epicyclic frequency for the Schwarzschild metric is 
given by k = ^. 
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Figure 10: Meridional cuts of the disk and the BH horizon in the xz-plane at several different phases f € (0 — 2n) of the radial 
oscillation. Each frame shows six density contours, equally spaced in logarithmic scale between the initial maximum density 
|Omax(0) and 10~ 6 p max (0). Also shown are the velocity field in the meridional plane and the location of the density maximum. 
The bottom two frames also show the position of the shock wave that is propagating outwards (thick black line). The last frame 
shows a brief episode of accretion. 
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of the disk bounces back earlier than the lower-density 
outer part. The collision of the former with the still- 
infalling lower-density outer material leads to formation 
of an outward-propagating shock wave (see two bottom 
panels on Fig. [10). The shock accelerates during prop- 
agation and reached relativistic velocities of ~ 0.3c in 
the rarefied outer disk shells. In order to quantify the 
amount of the shock dissipation, we evolve this model 
also with isentropic polytropic EOS with the same T — 
4/3 (in addition to the evolution with the non-isentropic 
T-law EOS), which does not allow entropy changes and 
thus shock heating. In this case, the oscillations exhibit 
little damping, while in the case of evolutions with the 
T-law EOS, the amplitude of |O max decreases by a fac- 
tor of ~ 2 in ~ 4 orbital periods. Note that dissipa- 
tion of oscillation kinetic energy into heat by shocks has 
been found to operate in many other astrophysical sce- 
narios, including damping of NS oscillations in a migra- 
tion from an unstable to a stable branch H1641 11651 , in 
phase-transition-induced collapse of NSs 11661 , as well 
as damping of ring-down oscillations of nascent proto- 
NSs formed in core-collapse supernovae (e.g., 1 167|) and 
accretion-induced collapse of white-dwarfs (e.g., 1 168 1). 



V. NON-AXISYMMETRIC INSTABILITIES 

In this section, we discuss non-axisymmetric insta- 
bilities in the disks. First, we test our method by re- 
producing results obtained by Kojima [80 1 for thick tori 
with negligible self-gravity on a Schwarzschild back- 
ground. In the next two sections we analyze the non- 
axisymmetric instabilities in our models first on fixed 
background, and then with fully dynamical general rel- 
ativistic treatment. 

In various astrophysical scenarios, the disk might 
be formed with different non-axisymmetric structures, 
leading to preferential excitation of specific unstable 
non-axisymmetric modes. In order to account for these 
possible scenarios, in our work we evolve initial disk 
models with added small non-axisymmetric perturba- 
tions at t — 0. 

We expect |63ll64"Il8"4l our thick disk models A and B 
to be dominated by non-axisymmetric unstable modes 
with azimuthal numbers m = 1 and m = 2, while 
for more slender model C instabilities with m = 3 and 
m = 4 should also play an important role. In our sim- 
ulations without initial perturbations, these modes are 
triggered by numerical errors and start to grow at a ran- 
dom moment in time. In order to explore evolution- 
ary scenario in which a particular mode is excited ini- 
tially, we add a small density perturbation of the form 
p = p[l + A cos m(q> — cpo)] to our initial disk models. 
We use a perturbation amplitude A — 0.001, which we 
have found to be both large enough to trigger the in- 
stability, and sufficiently small to remain in the linear 
regime. The amount of constraint violations due to these 
artificial perturbations, which is quickly suppressed by 
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Table III: Physical parameters of the initial disk models, used 
for comparison with calculations by Kojima 1 80 [ . Here, r c / r s is 
the location of the density maximum in units of the BH grav- 
itational radius r ? , r_ (r+) is the inner (outer) radius of the 
torus, and p max is the maximum density. For all models, the 
mass of the BH is M-q, specific angular momentum has con- 
stant value i = 4.0, and polytropic constant K = 0.06 (the lat- 
ter two quantities are given in the normalized system of units, 
in which G = c = Mq —!)• 



our constraint damping scheme, is too small to signifi- 
cantly affect the subsequent evolution of the disk. 

For each initial disk model A, B and C we have com- 
pleted a sequence of evolutions both in Cowling and 
in full GR treatment. As explained in the end of Sec- 
tion [nTC corresponding simulations are referred to by 
two- or three-letter notation (such as AC, AC1, AC2, 
AF, etc.), in which the first letter denotes the initial disk 
model (A, B or C), the second letter is either 'C for Cowl- 
ing or 'F' for full GR treatment, and the third one is the 
azimuthal number m of the initial perturbation (absent 
if evolved without initial perturbation). 



A. Comparison with previous work 

Below we present the growth rates of PP instability in 
evolutions of equilibrium tori in the Cowling approx- 
imation and compare them to the results obtained by 
Kojima llBTH , who studied PP instability in thick disks 
around Schwarzschild BHs using a linear perturbative 
approach. Kojima analyzed tori with constant distri- 
bution of specific angular momentum, constructed in 
a relativistic framework using the Abramowicz-Sikora- 
Jaroszynski (AJS) prescription |81|. Kojima calculated 
the m = 1 mode growth parameter 1/2 for several se- 
quences of disk models with different values of specific 
angular momentum I (= 3.8, 4.0, 4.2). In particular, 
Il80l found that, due to relativistic redshift effects, yi for 
GR models is generally smaller than in the Newtonian 
case 11631 . 

We have constructed a sequence of AJS tori of vary- 
ing radial extent with the same physical parameters as 
in [80 J. All these models have polytropic index T of 4/3 
and specific angular momentum t of 4.0. The mass of 
the BH is set to Mbh ~ 1- The models in this sequence 
are labelled as K1-K6, and their parameters are listed in 
Table lITl] To excite the m — 1 mode, we add a small 
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Model 
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Table IV: Parameters of the m = 1 PP instability, measured 
for the dynamical evolutions of Kojima disk models K1-K6. y\ 
and 1/2 are the pattern speed and growth rate parameters, Q p is 
the mode pattern speed, and r cr is the mode corotation radius. 



non-axisymmetric perturbation as described above, and 
evolve the disk in the Cowling approximation, using a 
curvilinear grid (as described in Section |IHC[ , adapted 
to each of the models. Parameters of the curvilinear 
grids for each model are listed in Table [TTJ We then mea 
sure the growth rate yi as described in Section HID 
above and compare it with the results of Kojima. 



B. Fixed background 

In this section, we analyze non-axisymmetric insta- 
bilities which develop when our initial disk models are 
evolved in Cowling approximation. We find that all of 
our models develop the PP instability. More specifically, 
for models A and B the fastest growing mode ism = 1, 
while for more slender model C it is m = 3. This is ex- 
pected from the Newtonian considerations |63 64, 84 1. 
Below, we first describe instabilities in models A and B, 
after which we focus on model C. As mentioned above, 
we evolve our models with and without artificial den- 
sity perturbations. For models A and B, we add m — 1,2 
and for model C we add m = 1,2,3 perturbations. No- 
tice that all simulations contain a spurious m = 4 per- 
turbation which is a numerical artifact of interpolation 
at the boundaries between four blocks near the equato- 
rial plane. The evolutions without artificial perturbation 
are therefore similar to the ones in which an m = 4 den- 
sity perturbation is added. 

To analyze the non-axisymmetric modes, we adopt an 
approach from Woodward et al [64J. Namely, we ex- 
pand the disk density in the e quator ial plane in a Fourier 
series, as explained in Section 
struct and analyze the following/owr diagrams: 
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Figure 11: Comparison of the m = 1 PP mode growth rates 1/2 
obtained by Kojima |80| in linear perturbative approach with 
the values measured from the evolutions of the same initial 
data models in Cowling approximation. The abscissa repre- 
sents the ratio of the inner radius of the disk r_ to the radius 
of the maximum disk density r c . Note that the error bars of 
the measured growth rates originate from the uncertainty in 
determining the time span of a clear exponential growth of the 
mode. 



Figure 1 1 shows the values of 1/2 as a function of the 
disk radialextent for our models and for those of Ko- 
jima. As we can see, 1/2 for our models are within esti- 
mated error bars from the values, calculated by Kojima 
in |80J, as it should be the case. 



(a) D m — t diagram shows logarithm of the normal- 
ized mode amplitude D m as a function of time at 
some radial location close to the disk density max- 
imum r c . The slope of this curve yields the growth 
rate 1/2- 

(b) D m — r diagram represents a radial profile of the 
normalized mode amplitude D m . This diagram 
will be helpful in identifying the type of non- 
axisymmetric instabilities [641. 

(c) cp m — t diagram displays a phase angle of the non- 
axisymmetric mode m as a function of time at a 
specified radius. The slope of this function deter- 
mines the mode pattern speed and parameter y\ 
that is related to it. 

(d) (p m — r diagram represents mode phase angle as 
a function of radius. This diagram also pro- 
vides a convenient way to identify the type of the 
mode 11641. In all our (p m — t diagrams, the disk 
rotates counterclockwise. 

We found that instabilities which develop during evo- 
lutions of models A and B are very similar. Therefore we 
present the properties of these instabilities on example 
of model A, while the case of model C will be described 
separately. 

Figure |i~2| shows the four mode diagrams for the case 
of model AC1, which represents time evolution of ini- 
tial disk model A in Cowling approximation with an 
added m = 1 non-axisymmetric density perturbation. 
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*For the simulation CF1, it was not possible to accurately determine the growth rate and pattern speed of the 
m = 1 mode. We classified this mode as a PP-type due to the character of its <p m — r and D,„ — r diagrams, 
which are typical for the PP-modes. 



Table V: Quantitative characteristics and types of the non-axisymmetric modes for the simulations studied in Section|V] For each 
of the simulations, the table lists one or two dominant unstable modes. 



Panel (a) shows the D m — t diagram for the first six non- 
axisymmetric modes with m = 1,2, ... ,6. Three lowest- 
m modes exhibit clear exponential growth, with m = 1 
being the dominant mode throughout the time span of 
the simulation (first ~ 6.5 orbital periods). This is the 
case because an m — 1 perturbation is artificially added 
from the beginning and has more time to grow and to re- 
main the dominant mode. The mode m — 2 has higher 
growth rate, but appears subdominant, since it is trig- 
gered later than the m = 1 mode and has less time to 
develop. It may eventually overshoot the m — 1 mode 
at a later time, when both modes reach nonlinear regime 
(not covered in this work, this will be a subject of our fu- 
ture publication). Notice that the relatively high values 
of the m = 4 mode amplitude are due to the effect of in- 
terpolation errors on four interblock boundaries of the 
grid near the equatorial plane. 



Panels (b) and (d) of Fig. 12 demonstrate the radial 
structure of the m = 1 mode at f = 4 t or b, when it is suffi- 
ciently developed. The D m — r diagram on Fig. 12 1 rep- 



resents the radial profile of the amplitude of the mode. 
This amplitude is highest near the edges of the disk, 
it does not have nodes (does not become zero), and 
it reaches its minimum near the radius of corotation. 
Previous works on PP instability in Newtonian grav- 
ity 1 61 . 84 J suggest that such radial behavior is a charac- 
teristic of the principal PP mode, or a mode of (0,0)-type 
in the classification of Blaes and Hawley [88J. The cp m — r 



diagram on Fig. 12 i has a specific S-shaped structure, 
which is also a well-known feature of the PP instability, 
discovered in previous Newtonian works [64, 84 J. 

Finally, panel (c) of Fig. 12 shows the cp m — t diagram 
for the m = 1 — 3 modes. iFshows that while all modes 
initially have arbitrary phases and pattern speeds, they 
eventually settle to the pattern speed of the dominant 
m — 1 mode possibly due to nonlinear interaction be- 
tween the modes. The pattern speed of the m — 1 mode 
is slightly below the speed of the disk at r c , which means 
that the mode corotation radius r cr lies just outside r c . 
The close proximity of the mode corotation radius to the 
radius of the disk density maximum is also typical for 
PP non-axisymmetric instabilities, as discovered in pre- 
vious Newtonian works lloTl 164 . 82 1 . All these features 
allow us to conclude that the observed m — 1 mode is 
indeed the PP instability. 



Figure 13 shows the set of four diagrams for simula- 
tion AC2, in which an m = 2 density perturbation is 
added initially. Panel (a) represents D m — t diagram for 
the modes with m = 1—6. In this case, only the mode 
with m = 2 exhibits pronounced exponential growth. 
The rest of the modes remain either stable or are not ex- 
cited, showing rapid growth only in the very end of the 
simulation, when the amplitude of m = 2 mode reaches 
nonlinear regime and coupling between the modes be- 
comes important. Panels (b) and (d) of Fig 13 show the 



radial profile and azimuthal shape of the mode. These 
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Figure 12: The behavior of the dominant non-axisymmetric modes and the radial character of the m = 1 (0,0) PP-mode in the 
simulation AC1: (a) D„, — t diagram shows the time evolution of the mode amplitudes D m for m = 1 — 6; (b) D m — r diagram 
shows the radial character of the m = 1 mode amplitude (in logarithmic scale); (c) <p m — t diagram shows the time evolution of the 
Fourier angle <p m , which allows to determine a pattern speed and a corotation radius for each mode; (d) <p m — r diagram shows 
the dependence of the mode Fourier angle from the radius r in equatorial plane at t/t orb ~ 4. The D m — r and (p m — r diagrams 
on panels (b) and (d) also show the locations of the mode corotation radius. 



are again typical for the principal (0,0)-type PP insta- 
bility with m — 2. Fig. 13 : shows the <p m — t diagram 
for the m = 1 — 3 modes. The dominant m — 2 mode 
rotates uniformly in the same direction with the disk, 
while other two modes do not exhibit clear rotation pat- 
tern until after ~ 4 orbital periods, when they align in 
phase with the m = 2 mode. The pattern speed of the 
dominant mode inferred from this plot corresponds to 
the corotation radius just outside of r c . Similarly to the 
simulation AC1 above, all these features are typical of 
the PP instability, studied in previous works in the New- 
tonian approximation [61 64, 82 J. 

While unstable modes in model B are very similar to 
those of model A, in model C higher-order PP modes 
become dominant. This is expected [63J since model C 
is more slender than models A and B. Indeed we ob- 
serve that the most unstable mode for the model C is 
the m = 3 PP mode, while m = 2 and m = 4 have com- 



parable but smaller growth rates. The four diagrams for 
m — 3, 4 modes are not qualitatively different from those 
for the m = 1,2 modes in simulations AC1 and AC2 
above. They also exhibit all the features typical of the 
PP instability described above. 

We can therefore conclude that the PP instability with 
various values of azimuthal number m is observed in 
all of our disk models in Cowling approximation. For 
reader's reference, the parameters of unstable modes 
calculated for all of our simulations are summarized in 
Tabled 



C. Dynamical background 

We now turn to the analysis of non-axisymmetric in- 
stabilities, which develop when the disks are evolved 
in a fully dynamical general relativistic framework. We 
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again consider evolutions with and without artificial 
density perturbation, adding m = 1,2 perturbations 
for models A and B and m — 1,2,3 perturbations for 
model C. While we observe PP type instabilities in the 
Cowling case, in the fully dynamical GR case we ob- 
serve two distinct types of instabilities: the Papaloizou- 
Pringle (PP) type and a GR analog of the so-called in- 
termediate type (I-type) instability |89, 91 J. Similarly to 
the Cowling case, instabilities in the moderately slender 
models A and B have very similar properties, so it suf- 
fices to present the result only for the case of model A. 
Model C is more slender and therefore favors instabili- 
ties with higher azimuthal numbers than those of mod- 
els A and B, so we consider this model separately. 

For the analysis of non-axisymmetric modes we adopt 
the same approach as in Section VB above for simu- 
lations on a fixed background, with one exception: for 
evaluating mode growth rates y%, instead of D m , we use 
quantities G m introduced in Section |IIID| above. This 



1' 

□de 



is necessary because the values of the mode amplitudes 
D m at a fixed radial location oscillate due to disk oscil- 
lations, making it hard to infer accurate growth rates of 



the modes. We have found that mode growth rate 1/2 
can be calculated more accurately from a time behavior 
of G m , because it is expressed in terms of integrals over 
the radius and as such it is much less affected by radial 
oscillations. With this exception, we follow the same ap- 
proach as described in Section |VB to determine mode 
types, growth rates and pattern speeds. These quanti- 
ties for all models are tabulated in Table IVl for reference. 

In most of the simulations with fully dynamical GR, 
the BH responds to the excitation of the m = 1 mode by 
developing an outspiraling motion, as described earlier 
in section [TV The position vector f bh of the BH starts 



to rotate with approximately constant angular velocity 
Ci^Hr while the length of this vector grows exponen- 
tially. In order to characterize this motion and study it 
in the context of the development of non-axisymmetric 
modes, we plot the time evolution of the BH position 
vector length y^h/ r g and phase angle cp^H on the G m — t 
and cp m — t diagrams, respectively. From these plots we 
can calculate the quantities an d yi (BH) that 

can be directly compared to those of non-axisymmetric 
modes. These quantities are also listed in the Table |Vj 
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Figure 14: The behavior of the dominant non-axisymmetric modes and the radial character of the m = 1 (0,0) PP-mode in the 
simulation AF1 (see captions to Fig. 12 for details). 



Figure [14] shows the time evolution and radial pro- 
files of the amplitudes and Fourier angles of the domi- 
nant unstable modes for the simulation AF1, which rep- 
resents fully dynamical GR evolution of the initial disk 
model A with an added m = 1 density perturbation. 
The top left panel contains the time evolution of G m for 
m = 1 — 4 and the normalized coordinate length of the 
BH position vector r^/rg. The diagram shows that the 
m — 1 mode is the dominant one. It also shows that the 
BH responds to the growth of the m = 1 mode and the 
distance from the BH to the origin grows exponentially 
at the same rate as the dominant m = 1 mode. 



Panel (c) of Fig. 



14 



shows the time evolution of the 
Fourier phase angles cp m for m = 1,2, 3, measured at a 
fixed radial coordinate location near the inner edge of 
the disk r_ . The phase of the dominant m = 1 mode 
after short initial readjustment exhibits almost uniform 
linear growth. Readjustment of the mode happens be- 
cause the added artificial perturbation initially does not 
have the right shape of the mode and needs some time 
(less than one orbital period) to readjust itself. Similar 



behavior is observed in the corresponding Cowling sim- 
ulation, but it is less pronounced (see the m = 1 phase 
angle during the first half orbit on Fig. 12 :). Phases of 



the rest of the modes initially do not show uniform lin- 
ear growth, but as the amplitude of the m = 1 mode 
increases, the higher order modes start to align in phase 
with the dominant m = 1 mode, most likely due to non- 
linear modes interaction. Notice that the pattern speed 
of the m = 1 mode calculated from this diagram is lower 
compared to the Cowling case. 



Panels (b) and (d) of Fig. 14 show the radial and angu- 
lar profiles of the dominant m = 1 mode at f = 2.5 t or jj. 
The location of r_, r+ and r c shown on the plot refers 
to the same time. Because the disk undergoes radial os- 
cillations, special care must be taken when calculating 
the mode corotation radius r cr , which is an important 
quantity that characterizes non-axisymmetric modes. To 
find T cr , we solve an equation between the mode pattern 
speed Clp and the disk angular velocity: Ci(r cr ) = Ci p . 
The latter changes during the evolution of the disk, so 
strictly speaking the value of r cr will also depend on 
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time. However, we have found that in all our simula- 
tions the change of the profile of O in the course of the 
evolution is very small, therefore within the measured 
accuracy the corotation radius is independent of time. 
The values of r cr / r c , calculated in this way for all sim- 
ulations on a dynamical background are also listed in 
Table [V] For simulation AF1, r cr lies outside of r c with 
J"cr/J"c ~ 1.17. Such a value of r cr / r c is typical for a PP in- 
stability and comparable to the values observed in pre- 
vious Newtonian studies with a moving central object 
and a massive self -gravitating disk 164] . 

Similar to the Cowling case, the radial profiles of D m 
and cp m show structural features that are typical for a PP 
instability of (0,0)-type in the Blaes and Hawley classifi- 
cation 1 88 1. In particular, the mode amplitude displayed 
on Fig. 14 r> is higher near the edges of the disk and has 
a minimum near corotation r cr . The mode amplitude 
does not vanish anywhere in the disk, which character- 
izes the mode as having type (0,0). The profile of the 
Fourier angle of the mode, shown on Fig. 14 i, has a spe- 
cific S-shaped structure, that consists of a trailing spiral 
pattern outside the corotation radius 8 , a leading spiral 
pattern inside the corotation radius and a short segment 
near r cr that connects the two spiral patterns. 

However, compared to the Cowling approximation, 
the growth rate of the m = 1 mode is amplified by al- 
most a factor of ~ 1.5, and there are reasons to believe 
that the outspiraling motion of the BH is responsible for 
this amplification. The mechanism which drives the un- 
stable outspiraling motion of the BH is similar to the one 
described in 11941 . The disk creates a hilltop potential, 
which has a maximum at the origin, so the BH initially 
is located at the point of unstable equilibrium. The BH 
can reduce its potential energy by converting it into ki- 
netic energy of orbital motion around the common cen- 
ter of mass (CM) of the disk+BH system. Such orbital 
motion requires angular momentum which can be bor- 
rowed from the disk through the development of a non- 
axisymmetric m = 1 mode. Because the orbital motion 
of the BH requires a compensating displacement of the 
disk CM, removing angular momentum from the disk 
increases the amplitude of the m = 1 mode, which in 
this case is the PP-mode. 

Next we consider the dominant non-axisymmetric 
modes that develop in the simulation AF2, in which 
an m = 2 density perturbation was added. Figure 15 
presents the time evolution and radial profiles of G m ancT 
q> m of these modes. Comparing Fig. 15 to Fig. 12 which 
presents the same set of diagrams for the CowTmg sim- 
ulation AC2, we can see that the type of the m = 2 mode 
in this simulation is quite different from the PP one ob- 
served in the simulation AC2. First, on the D m — r dia- 



gram in Fig. 15 :>, the minimum of D m lies inside r c , un- 



like the case of the PP mode where such minimum is 
located close to the mode corotation radius r cr . Second, 
the (p m — r diagram of the m = 2 mode in Fig. 15 1 con- 
sists only of trailing spirals and does not have aleading 
spiral pattern inside corotation, as it would be the case 
for PP modes. Finally, as Fig. [l5p l illustrates, the mode 
phase angle makes an abrupt turn by n/ 2 radians near 
r c , which would be the case when the disk were sub- 
jected to an elliptic (bar-like) deformation. 

Overall, this mode looks very similar to the so-called 
intermediate type (I-type) modes found in earlier stud- 
ies using Newtonian gravity (see [64, 89-91]). In partic- 
ular, 1 64 1 observed the I-modes in their 3D Newtonian 
simulations of self-gravitating disks with various disk- 
to-central object mass ratios Mq / M c and values of pa- 
rameter T/|W| (see Section 4 in [64J). A subset of their 
models that develop the I-mode instability (namely E31 
and E32) have parameters M D /M C = 0.2, T/\W\ ~ 0.47 
and ArcGp/dc ~ 4, which are comparable to those of 
our models A or B (listed in Table |IJ. Notice that in this 
simulation the m = 1 mode is also excited, as can be in- 



ferred from Fig. 15 i. The growth rate and pattern speed 
of this mode is me same as in simulation AF1 (cf. Ta- 
ble]^. The m = 1 mode development is again correlated 
in time with o utsp iraling motion of the BH, which is ap- 
parent in Fig. |15fr . The analysis of the mode character, 
similar to the one performed above confirms that in this 
simulation the m = 1 mode has the same (0,0) PP type as 
in the simulation AF1. It grows faster than m — 2 mode, 
so that both modes become of comparable amplitude by 
the end of the simulation. The radial character of the 



We remind the reader that the disk rotates counterclockwise on all 
cp m — r diagrams. 



two co-existing modes is also depicted on Fig. 16 which 
shows a sequence of snapshots of the <p m — r diagrams 
at different times, superimposed with the disk density 
in equatorial plane. 

In the case of model C, we used artificial density per- 
turbations with m — 1, 2, 3 and observed the develop- 
ment of four unstable modes. The same analysis that 
we have done before for models A and B reveals that 
only the m = 1 mode has PP type, while the modes 
with m = 2, 3, 4 have I-type. The I-modes with m = 3 
and m = 4 represent triangular and square deforma- 
tions of the disk. Such modes were previously observed 
in Newtonian simulations of narrow self -gravitating an- 
nuli Il90l . Figure 17 shows the time behavior of G m and 
r BH I r g for simulations of model C with different added 
perturbations. As can be seen from these plots, the 
fastest growing mode is the I-mode with m = 3, while 
the slowest one is the m = 1 PP mode (cf. Table |V) , 
The latter remains subdominant in all simulations, its 
growth does not show clear exponential behavior and it 
is poorly correlated with the motion of the BH. This hap- 
pens because the m — 1 mode does not form a global 
coherent pattern. In this case, the quantity G m does 
not correspond to the amplitude of a global mode, but 
rather to some combination of local m = 1 Fourier har- 
monics. However, as soon as the amplitude of the m = 1 
mode reaches G m ~ 10~ 4 , it shows a better correlation 
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Figure 15: The behavior of the dominant non-axisymmetric modes and the radial character of the m = 2 1-mode in the simulation 
AF2 (see captions to Fig. 12 for details). 



with the motion of the BH, as can be seen on Fig. 17 for 
simulations CC1 and CC3. 

In model C, the m — 1 mode starts growing when the 
other modes with higher m have already reached sig- 
nificant amplitudes (G m ~ 0.1). Moreover, the growth 
rate of the m = 1 mode has different values depend- 
ing on which of the higher-m modes dominates the dy- 
namics of the disk. For example, in models CF1, CF2, 
and CF3, in which the modes with m = 4, m = 2 and 
m — 3 reach the largest amplitudes, the m — 1 mode 
has the growth rates of 0.11, 0.23 and 0.28, respectively. 
Such behavior of the growth rates is likely to be a result 
of non-linear interaction between these modes. In other 
words, m = 1 mode is never linearly independent from 
the m > 1 modes, making it hard to explore unambigu- 
ously the influence of the BH motion on the growth rate 
of the m = 1 mode. 

Finally, we need to point out that the comparisons 
with the Cowling simulations presented in this section 
should be taken with a grain of salt. The initial per- 



turbation introduced by the blending of two metrics in 
the initial data changes the BH mass by a small amount 
(w 2.5%), which not only drives the disk out of equi- 
librium, but also changes the equilibrium configuration 
itself. Therefore, strictly speaking, the disks in Cowling 
and full GR represent different objects and cannot be di- 
rectly compared. As a result, these disks will have dif- 
ferent evolutionary paths not only due to dynamical GR 
effects, but also due to differences in BH masses. How- 
ever, because the latter is small, we believe that the dif- 
ferences in evolution are mainly caused by the former, 
while the latter should not affect the time evolution sig- 
nificantly. For example, for models A and B which have 
different disk-to-BH mass ratios (0.24 vs 0.17), we do not 
observe qualitative differences in time evolution, and 
quantitative differences are small (e.g. the differences in 
the growth rates are within 6%). Therefore, we conclude 
that most of the differences between our Cowing and 
full GR simulations are caused by the effects of GR. 
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Figure 16: A sequence of four successive snapshots of the disk density in the equatorial plane, combined with the corresponding 
cp m — r diagrams for the m = 1 and m = 2 modes for the simulation AF2. The radial character of the m = 1 and m = 2 modes is 
represented by a sequence of red and blue dots resp. 
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Figure 17: G m — t diagrams for modes with m = 1 — 4 in simulations with disk model C 
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Source 


LIGO 


Adv. LIGO 


ET 


A at 10 kpc 


1 


1 


1 


C at 10 kpc 


2 


1 


1 


A at 18 Mpc 


1.7 x 10 6 


7900 


40 


C at 18 Mpc 


6.4 x 10 6 


80000 


150 



Table VI: Minimal number of wave cycles needed for grav- 
itational waves from non-axisymmetric instabilities to be de- 
tectable by LIGO, Advanced LIGO and Einstein Telescope 
(ET). It is assumed that the value of the amplitude of the in- 
stability is at least D m = 0.1 during this time. Estimates are 
given for the sources with parameters of models A and C, lo- 
cated at distances 10 kpc and 18 Mpc. 



D. Gravitational Wave Detectability 

As mentioned before in Section |IV| all of our disk 
models are unstable to non-axisymmetric modes. Once 
formed, these modes start growing exponentially until 
they reach a saturation regime due to non-liner effects. 
This process is accompanied by a redistribution of the 
angular momentum of the disk until the profile of the 
specific angular momentum becomes steep enough for 
the disk to be stable to non-axisymmetric instabilities 
(see related discussion in [86 91 J). Before the angular 
momentum is redistributed, and even after the disk be- 
comes stable, the amplitude of non-axisymmetric modes 
in the disk is likely to remain high (possibly near the sat- 
uration level) 1 86]. The presence of non-axisymmetric 
deformations in the disk leads to emission of potentially 
detectable gravitational radiation. 

Below we give estimates of the detectability of the GW 
signal from saturated non-axisymmetric instabilities in 
our disk models. We make our estimates based on the 
Newtonian quadrupole formula for the initial disk mod- 
els with added m — 2 mode with an amplitude D m — 
0.1. We calculate an approximate number of cycles that 
the instability needs to remain at that amplitude in or- 
der for the emitted GW to be detectable. These numbers 
are listed in Table |Vlj for an event at a distance of 10 kpc 
(our galaxy) and 18 Mpc (a distance to the Virgo cluster), 
for models A and C. The table shows that an event in 
our galaxy will be detectable with the current LIGO de- 
tector even with a single cycle of the non-axisymmetric 
mode. An event in the Virgo cluster, on the other hand, 
is unlikely to be detectable with the current LIGO detec- 
tor, since it would require unrealistically large number 
(> 10 6 ) of cycles. Second and third generation detectors 
such as the advanced LIGO and the Einstein Telescope 
can detect such events if non-axisymmetric modes per- 
sist for ~ 10 4 — 10 5 and ~ 40 — 150 cycles, respectively. 
Finally, we point out that it is currently unclear how 
long a non-axisymmetric mode in a given disk model 
will persist in non-linear regime. This is likely to depend 
on the details of non-linear mode properties, accretion 



rate, magnetic fields and the thermodynamic state of the 
disk matter Il86l . 



VI. CONCLUSION 

In this paper we have explored non-axisymmetric in- 
stabilities in self-gravitating disks around black holes 
(BHs) using three-dimensional hydrodynamical simula- 
tions in full general relativity (GR). We studied several 
moderately slender and slender models with disk-to-BH 
mass ratio ranging from 0.11 to 0.24. The parameters of 
these models are listed in Table |l] 

To obtain a self-consistent equilibrium disk model 
outside BH, we solve the coupled system of Einstein 
constraints and Euler equations using an iterative Green 
functions approach, implemented in the RNS code [150 1. 
To avoid coordinate singularities, we transform the sta- 
tionary initial data outside the BH horizon from quasi- 
isotropic to non-singular horizon-penetrating coordi- 
nates. We set the data inside the BH horizon to the an- 
alytic Kerr-Schild solution and smoothly blend it with 
the computed data outside the horizon. 

We evolve the metric using a first-order form of the 
generalized harmonic formulation of the Einstein equa- 
tions with adaptive constraint damping. The metric 
evolution equations are discretized on multiblock grids 
and solved using 8th order finite difference operators. 
We evolve the matter with relativistic hydrodynamics 
equations in flux-conservative form, using a finite vol- 
ume cell-centered discretization scheme. We use a T-law 
equation of state to model disk matter. Our numerical 
approach makes extensive use of the curvilinear mesh 
adaptation in order to achieve desired resolutions in dif- 
ferent parts of the domain. 

We did not observe the runaway instability in our 
models, which could have developed due to the disk 
overfilling its toroidal Roche lobe in the process of radial 
oscillations. Such radial axisymmetric oscillations of the 
disk around its equilibrium state are triggered in all of 
our simulations by an axisymmetric perturbation in the 
metric due to the blending of two metrics inside and 
outside the BH horizon (see Section III A I. Although our 
initial disk models are close to overfilling the toroidal 
Roche lobe, the radial oscillations do not lead to the de- 
velopment of the runaway instability within several ini- 
tial orbital periods that we have simulated. However, 
this result may be specific to the particular model that 
we focused on; we can not exclude the possibility that 
the runaway instability develops in models with differ- 
ent initial parameters. 

In all models that were explored we observed unsta- 
ble non-axisymmetric modes. We have performed de- 
tailed analysis of these modes to determine their types, 
growth rates, radial profiles and pattern speeds (see Ta- 
ble|V|. For all simulations in the Cowling approximation 
we observe the development of the Papaloizou-Pringle 
(PP) instability with m = 1 — 4. In this case, the az- 
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imuthal number m of the fastest growing mode depends 
on the disk slenderness: for moderately slender models 
A and B such mode is m = 2, while for more slender 
model C, it is m = 3. In the simulations in full GR, we 
observe two distinct types of instabilities. The unstable 
mode with m = 1 has PP type, similar to the one ob- 
served in Cowling case. Unstable modes with m > 1 
become the intermediate modes (I-modes), representing 
elliptical, triangular or square deformations of the disk. 
In full GR, the fastest growing mode is m — 1 in models 
A and B, and m = 3 in model C. 

In the full GR case, the development of the m = 1 PP 
mode is accompanied by an outspiraling motion of the 
BH. The distance from the BH center to its initial posi- 
tion has the same growth as that of the m = 1 mode am- 
plitude. We find that due to this motion, the growth rate 
of the m = 1 mode is amplified by a factor of « 1.5 com- 
pared to the Cowling case for massive models A and B. 
This amplification makes the m — 1 PP mode the fastest 
growing one in the models A and B, while in the case of 
less massive model C, this mechanism is not as efficient. 
The overall picture of the unstable modes in full GR is 
qualitatively similar to and consistent with the Newto- 
nian case 1611 1801 1901. 

Evolution of non-axisymmetric instabilities in non- 
linear regime will be associated with the emission of 
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high-frequency gravitational radiation. In Table 
we give rough estimates of the detectability of this ra 
diation in terms of the number of cycles that a non 
axisymmetric deformation must persist in order to be 
detectable. While even a single cycle of gravitational ra- 
diation from this deformation is detectable if occurs in 
our galaxy, for more reasonable distances such as Virgo 
cluster it is only detectable with Advanced LIGO, and 
only in the case that the disk deformation persists for 
thousands of cycles (see Table |Vj). It is currently unclear 
how long a non-axisymmetric deformation can persist 
in non-linear regime. 

Finally, we would like to point out limitations of our 
current simulations. We use simplified initial disk mod- 
els and do not include realistic microphysics, neutrino 
cooling and magnetic fields. Future studies of non- 
axisymmetric instabilities should take into account these 
effects, as well as consider larger set of parameters, such 
as non-constant angular momentum distribution, vari- 
ous disk sizes, masses, BH spins etc. The properties of 
the disk in the nonlinear regime, such as the persistence 
of non-axisymmetric structures in realistic disk models 
should also be addressed. 



VII. ACKNOWLEDGEMENTS 

The authors would like to thank our colleagues Eloisa 
Bentivegna, Peter Diener, Juhan Frank, Tyler Landis, 
Frank Loffler, Luis Lehner, Christian D. Ott, Jorge 
Pullin, Jian Tao, Manuel Tiglio, and Joel Tohline for 
valuable discussions and ideas. We also thank Ya- 



sufumi Kojima for providing his data for comparison 
with our results (Fig. [TT) . This work is supported 
by the NSF grants 072M5 (Alpaca), 0904015 (CIGR), 
and 0905046/0941653 (PetaCactus). The simulations 
were performed using the supercomputing resources 
Ranger and Lonestar at TACC via the NSF TeraGrid, and 
Queenbee at LONI. We also used the PetaShare infras- 
tructure to store the data from our simulations. N. S. 
acknowledges the hospitality of the University of Tubin- 
gen, and O.K. acknowledges the hospitality of AEI nu- 
merical relativity group. 



Appendix 

1. Transforming initial data to horizon-penetrating 
coordinates 



Here we present the transformation of the station- 
ary axisymmetric initial data in quasi-isotropic coordi- 
nates to time-independent horizon-penetrating coordi- 
nates which is used in Section [Till m order to remove 
the degeneracy at the BH horizon. In our initial data, 
the metric is given in general form |[6| and represents an 
axisymmetric deformation of a Schwarzschild BH by a 
massive equilibrium torus. The sought transformation 
has to satisfy the following requirements: 

• the metric in the new coordinates is time- 
independent; 

• the metric does not have pathologies (degeneracy 
or divergence) at the event horizon; 

• the three-metric on f = const, foliation is positive 
definite (i.e. the t = const, foliation is spacelike). 

We build our transformation by analogy with the 
transformation from isotropic Boyer-Lindquist coordi- 
nates 1 169 1 to horizon-penetrating Kerr-Schild coordi- 
nates [170 J of Schwarzschild spacetime (see also 111 71 II 
for a general case of rotating BH). In this case, the line 
element has the following form: 



dsf 



1 - m/2r* 
l + m/2r* 



dt 2 + rp i (drl+rldCi 2 ) 



where m is the BH mass, r* is the isotropic radius, dCi 2 = 
d8 2 + sin 8 2 dq> 2 is the solid angle element, and ip = 1 + 
27- is the conformal factor. At the event horizon ; 7 = 
m/2 the determinant of the isotropic metric is zero. 

In the horizon-penetrating Kerr-Schild coordinates, 
the line element will be: 



dsf s = - (1 - H)dt z + 2Hdtdr + (1 + H)dr z + r 2 dfi 2 

where H = 2m/r and r = r*(l+m/2r*) is the 
Schwarzschild radial coordinate. 
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The Jacobian of the transformation from the isotropic 
(t,r*,6,(p) to the horizon-penetrating coordinates 
(t, r, 9, cp) has the following form: 



D(t, r*,9,cp) 



1 


H 
1-H 











!'* 








ry/l-H 








1 














1 



D{t,r,9,q>) 



The metric in the new coordinates remains indepen- 
dent of time, which is also true for an arbitrary transfor- 
mation with the following Jacobian: 



D{t,r„9,q>) 
D(t,r,9,q>) 



1 f(r,9) h{r,6) 

g{r,6) p(r,8) 

10 

1 



(Al) 



where the functions {/, g,h,p} do not explicitly depend 
on time and are only constrained by the regular Jacobian 
integrability conditions: 



df(r,6) _ dh(r,9) 
39 ~ dr ' 



dg(r,9) = dp(r,9) 
d9 dr 



For the case of a general axisymmetric spacetime, the 
metric in quasi-isotropic coordinates is given by: 



gtt 

e 2a 



-vg<p 9 o 



r 2 e 2a 



-vg<p<p 





where gtt = —A 2 



and gipm = B 2 A 2 r 2 sin# 2 . 



We can select the set of functions \f,g,h,p} in the Ja- 
cobian (Al) above to construct a transformation from 
the quasi-isotropic to horizon-penetrating coordinates, 
which satisfies the above requirements for an arbitrary 
stationary axisymmetric deformation of Schwarzschild 
if we choose: 

f(r,6) = 1 - 1/A 2 (r,0), g(r,6) = e~ a ^ / 'A(r,0). 

Since the metric potentials A(r*,6) and a(r*,#) are 
known as functions of r* and not r, we need to express 
the new radial coordinate r in terms of r* . The required 
relations in the following differential form are obtained 
by inverting the Jacobian (Al): 



dr* 
dr 



?=const. 



dr* 
~d~6 



0. 



r*=const. 



These need to be integrated along the radial coordinate 
from r* /, to r* for each 9: 



r(r*,0) - r h 



where r;, is the radius of the event horizon in new coor- 
dinates. The remaining two unknown functions h and p 
can be calculated by ID integration of the Jacobian inte- 
grability conditions: 



h(r,6) 
p(r,0) 



C dr' 

Jr h 


df(r',9) 


d9 


f dr' 

Jr h 


dg(r',8) 
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-J" 





r. \' e {UV m 



di. 



After the transformation, the metric in the new 
horizon-penetrating coordinates has the following form: 



gks 



gtt 



e 2a g 2 



fgtt hgtt 

gttf fhgtt + e 2a gp -ufg v<P 
h 2 gtt + e 2a (p 2 + rl) - cvhg f(p 



where ellipsis indicate matrix elements which can be 
filled in by symmetry. In the limit of r — >■ the func- 
tions / and g tend to infinity, while h and p vanish. The 
resulting metric at the horizon remains finite and non- 
degenerate: 



lim gks 

T-tTf, 





C 



C e 2a r 2 





where C = lim,-^^ [fhgtt + s 2a gp] is a finite constant. 

2. Stable evolution of a uniformly rotating polytrope 

In this appendix, we present results of testing the time 
evolution of a uniformly rotating polytropic star for nu- 
merical stability and convergence. In our tests, we use 
geometrized units based on the solar mass, in which 
G = c = .M q = 1. The parameters of the star in the ge- 
ometrized and CGS units are summarized in Table I VII I 
We use a thirteen-block cubed sphere system that was 
described in 1 102 1 (see Fig. [lb and the related discussion 
in Section II A above). For the current setup, we fix the 
sizes of the blocks by choosing r^ = 2.5, ri = 9 and 
f2 = 14 (see Section 4 of [102] for definition of r^, r\ and 
ri), with each block having an equal number of N 3 grid 
cells. The sizes of the domain and its blocks are selected 
in such a way that the star occupies ~ 90% of the en- 
tire domain in radial equatorial direction, and the inner 
seven blocks of the system lie inside the star. This setup 
allows to test how much the accuracy and convergence 
of our numerical scheme are affected by interpolation er- 
rors on the interblock boundaries which thread the bulk 
of the star. 
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geom. 


CGS 


polytropic scale K 


100 


1.46 • 10 5 cm 5 g" 1 s" 2 


polytropic index T 


2 


2 


central rest-mass density p c 


0.001 


6.17 ■ 10 14 g cm" 3 


ratio Rp/R e 


0.7 


0.7 


ADM mass M 


1.49 


1.49 


rest mass Mg 


1.59 


1.59 


equatorial radius R e 


12.32 


1.823 ■ 10 6 cm 


angular momentum / 


1.32 


1.16- 10 49 gcmV 2 


normalized ang. mom. // M 2 


0.59 


0.59 


kinetic / binding en. T/\W\ 


0.0748 


0.0748 


angular velocity Q 


0.0215 


4300 s" 1 


Keplerian angular velocity 


0.0286 


5801 s" 1 


rotational period P 


292.1 


1.44 ■ 10~ 3 s 



Table VII: Physical parameters of the uniformly rotating poly- 
tropic star used for the code tests, in geometrized and CGS 
units, where: Rp / R e is the ratio of the polar to equatorial radii 
of the star, // M 2 is its angular momentum, normalized with 
the square of the ADM mass of the star M, and T/\W\ is the 
ratio of the kinetic to binding energy of the star. 



The stability of the numerical scheme for evolving 
the spacetime metric depends on the numerical dissipa- 
tion parameter e [144J and the constraint damping co- 
efficients k, 72 (see Section [IIE[ . In general, higher val- 
ues of numerical dissipation restrict the timestep, while 
lower values are undesirable because they do not pro- 
vide enough suppression of the numerical noise, which 
needs to be dissipated for stability 11721 . For the current 
setup, we choose e — 0.2 and k = 72 = 0.1. Values of the 
constraint damping parameters higher than ~ 0.5 lead 
to numerical instabilities in our simulations of stars. 

Initial data for the time evolution is generated by the 
RNS code fT50| , which uses the KEH(SF) method 11501 



11521 to produce equilibrium models of stationary rotat- 
ing relativistic stars. Since RNS is a 2D solver which 
uses its own grid that is different from the 3D multi- 
block grid of our time evolution code, we interpolate the 
data from 2D grid to the 3D multiblock grid using 4-th 
order Lagrange interpolation. Also, because the vari- 
ables in the GH formulation contain first derivatives of 
the metric and because the resolution on the 2D grid is 
usually much higher than on the multiblock 3D grid, we 
perform numerical differentiation on the 2D grid. The 
resulting derivatives are then interpolated onto the 3D 
grid. Note that the interpolation procedure is not con- 
sistent with the Einstein constraint equations, and hence 
produces numerical noise. 

The system is evolved up to f = 350, which cor- 
responds to 1.73 ms, or 10 dynamical timescales of 
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Figure 18: Time evolution of Lj norms of the Hamiltonian 
constraint (top panel) and density solution error Sp/p max (0) 
(bottom panel) for three different resolutions. 



the star 9 . In vacuum regions outside the star, we use 
an artificial atmosphere, which has density of jO a tm = 
10 _7 l o max (0), where p m ax(0) is the maximum density at 
t = 0. If during the evolution the density in a cell drops 
down below a threshold value set to p^ = 2 p a t m , the 
density in this cell is reset to the artificial atmospheric 
value. To estimate the accuracy of our code, we have 
performed a convergence study using three different 
resolutions with N 3 = 20 x 20 x 20, 40 x 40 x 40 and 
80 x 80 x 80 grid points in each block. We have ana- 
lyzed various integral norms of the errors in all evolved 
variables, including 50 spacetime variables, 5 primitive 
variables and 5 conserved variables. We have also an- 
alyzed integral norms of the Hamiltonian and momen- 
tum constraints, as well as the behavior of conserved in- 
tegral quantities such as total rest mass and total angular 



9 The dynamical time to is defined as to = R e v 7 R e I M, where R e is 
a proper equatorial circumferential radius, and M the ADM mass 
of the star. It corresponds to the inverse of the orbital frequency 
fl = a/M/R! at R e . 
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Figure 19: Time evolution of the total rest mass (top panel) and 
the total angular momentum (bottom panel) for three different 
resolutions. Both quantities are normalized to their values at 
t = 0. 



momentum. 

In all cases, we observe the expected 2-nd order con- 
vergence. As an example, Fig. 18 (top panel) shows the 
time evolution of the L\ norm of the normalized den- 
sity deviation Sp = [p(t) — p(0)] /p m ax{0) for the three 
resolutions. Due to accumulation of truncation errors, 
this deviation exhibits a steady growth (modulo small 
variations because of oscillations of the star) throughout 
entire evolution. The deviation for N = 40 is larger than 
that for N = 80 by a factor of « 4, which is a clear signa- 
ture of 2-nd order convergence. However, the deviation 
for N = 20 is larger than that for N = 40 by a smaller 
factor of « 1.5, which means that the resolution N = 20 
is insufficient for achieving a convergent regime. A sim- 
ilar convergent behavior is observed for integral norms 
of the deviations of all of the rest of the variables. 

Figure [18] (bottom panel) shows the plot of the L\ 
norm of trie Hamiltonian constraint violation as a func- 
tion of time. This quantity is not zero at f = 0, since 
initial conditions were interpolated from the 2D grid 
and interpolation errors were introduced. However, be- 
cause of the constraint damping scheme, the Hamilto- 
nian constraint violation significantly drops for medium 



and high resolutions within the first 0.2 ms. During 
subsequent evolution the value of the Hamiltonian con- 
straint remains stable and clearly shows 2-nd order con- 
vergence with resolution, i.e. the values of the Hamilto- 
nian constraint for N = 20, 40 and 80 are in an approx- 
imate ratio 16 : 4 : 1. Momentum constraints show a 
similar behavior. 



Figure 19 demonstrates time evolution of the total rest 
mass (upper panel) and total angular momentum (lower 
panel) of the star. In our numerical simulations these 
quantities are not conserved mostly due to interpolation 
errors on the interblock boundaries that pass through 
the bulk of the star. By the end of the simulation, for 
N = 20, 40 and 80, the total rest mass decreases by 
0.88, 0.84 and 0.25 percent, while the total angular mo- 
mentum decreases by 0.07, 0.12 and 0.05 percent. These 
numbers show that the smallest necessary resolution for 
the convergent regime is N = 40, which amounts to 
» 70 — 100 points across the star. 



3. Fundamental modes of a TOV star 

As another test of the coupling between the GR 
and hydro parts of the code, we evolved a Tolman- 
Oppenheimer-Volkoff (TOV) solution on a seven-block 
system, and measured the frequencies of its fundamen- 
tal oscillations both in the Cowling approximation and 
in full GR. In these tests, we use geometrized units in 
which G — c — Mq — 1. We choose a star with T = 2, 
K — 100 and the value of rest-mass density in the center 
p c = 1.28 ■ 10~ 3 . These parameters produce a TOV star 
with gravitational mass M = 1.4 and circumferential ra- 
dius R e = 9.8. This system has already been extensively 
studied in the literature and used for the assessment 
of relativistic hydrodynamical codes (e.g. B1641 11731 ). 
The seven-block cubed sphere system that we used (see 
Fig. [l^) has the outer radius R = 12, which makes the 
star occupy 82% of the domain in radial direction and 
leaves extra room for small oscillations. The size of the 
cubical block in the center is a — 4.8, placing it com- 
pletely inside the star. The bulk of the star is there- 
fore threaded by interpolation boundaries between the 
blocks. The cubical block contains N 3 volume cells, and 
the outer blocks have N 2 x (2N) cells. For the tests, we 
used resolutions N = 20, 40 and 80, which roughly cor- 
respond to 40, 80 and 160 points across the star. 

To observe and measure the fundamental mode, we 
artificially add a small initial perturbation, roughly cor- 
responding to the shape of the mode: 

Sp , nr 
-f = Acqs— , 
P 2R e 



The amplitude was chosen to be A = 0.005. Fig. 20 
displays the resulting oscillatory behavior of the rest- 
mass density in the center of the star for three different 
resolutions. Left and right panels correspond to fixed 
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Figure 20: Time evolution of the density p c (t) at the center of a TOV star, normalized by its initial value po, for three different 
resolutions. Left panel: Cowling approximation. Right panel: fully dynamical GR case. 
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Figure 21: Power spectrum of the density oscillations p c {t) 
at the center of a TOV star for the highest resolution simula- 
tions in Cowling approximation (red solid line) and in full GR 
(green solid line). Also shown are derivatives of the power 
spectrum with respect to the frequency obtained using the 
central finite-differencing scheme. The derivatives allow to lo- 
calize peaks in the power spectrum more accurately. Vertical 
axis has arbitrary units. 



(Cowling approximation) and dynamical spacetime ge- 
ometries, respectively. Oscillations of the density are ac- 
companied by a secular drift, which reflects accumula- 
tion of truncation errors and converges away with reso- 
lution at approximately second-order convergence rate. 
Conserved quantities such as the total rest mass and 
the total angular momentum (not shown) also exhibit 
the second-order convergence, as expected. In particu- 
lar, for the Cowling approximation case, the simulation 
continued up to 20 ms, and the final rest mass is con- 
served up to 8.3, 3.7 and 1.1 percent for resolutions with 
N = 20, 40 and 80. For the full GR case, the simula- 
tion continued for 6 ms and the rest mass is conserved 



up to 3.2, 1.3, 0.4 percent for the same three resolutions. 
This shows that the rate of the mass loss in Cowling and 
full GR simulations is approximately the same, as ex- 
pected. Because the bulk of the star is threaded by inter- 
polation boundaries between the blocks, the mass loss is 
significantly higher in this setup than in case of a regu- 
lar Cartesian grid, where we normally observe that the 
mass is conserved up to 7-8 significant digits for a simi- 
lar resolution. 

A Fourier transform of p c {t) allows to measure the 



Fig. 
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frequencies of the dominant oscillation modes 
shows the Fourier power spectrum of Pc(t) in Cowl 
ing and full GR cases for simulations with the high- 
est resolution N = 80. Both spectra contain three eas- 
ily identifiable peaks corresponding to the fundamen- 
tal radial modes F, H\ and Hz- The same plot also 
shows derivatives of the spectral power with respect to 
the frequency, computed using the central finite differ- 
encing scheme. Zeroes of these numerical derivatives 
provide accurate estimates of the location of frequency 
peaks. The frequency of the F-mode in Cowling approx- 
imation is v(F) = 2.684(40) kHz, which is in agree- 
ment with the value of 2.706 kHz, found in [164]. In the 
fully general relativistic case, we obtain the frequency 
v(F) = 1.440(50) kHz, which also agrees with the value 
1.458 kHz, found in Ifl63| . Note that the error in the 
values of fundamental frequencies above is estimated 
as the distance between the root of the power spectrum 
derivative and the nearest point with a non-zero value. 



31 



[1] S. E. Woosley, ApJ 405, 273 (1993). 

[2] S. E. Woosley and J. S. Bloom, Annual Rev. Astron. As- 
trophys. 44, 507 (2006), astro-ph/0609142. 

[3] Y. Sekiguchi and M. Shibata, ArXiv e-prints (2010), 
1009.5303. 

[4] D. Proga, A. I. MacFadyen, P. J. Armitage, and M. C. 

Begelman, ApJ 599, L5 (2003), arXiv:astro-ph/0310002. 
[5] S. Fujimoto, K. Kotake, S. Yamada, M. Hashimoto, and 

K. Sato, ApJ 644, 1040 (2006), arXiv:astro-ph/0602457. 
[6] L. Dessart, A. Burrows, E. Livne, and C. D. Ott, ApJ 673, 

L43 (2008), 0710.5789. 
[7] E. O'Connor and C. D. Ott, ArXiv e-prints (2010), 

1010.5550. 

[8] S. Nagataki, ApJ 704, 937 (2009), 0902.1908. 

[9] S. Harikae, T. Takiwaki, and K. Kotake, ApJ 704, 354 

(2009) , 0905.2006. 

[10] D. Lopez-Camara, W. H. Lee, and E. Ramirez-Ruiz, ApJ 

692, 804 (2009), 0808.0462. 
[11] A. I. MacFadyen and S. E. Woosley, ApJ 524, 262 (1999), 

arXiv:astro-ph/9810274. 
[12] A. I. MacFadyen, S. E. Woosley, and A. Heger, ApJ 550, 

410 (2001), arXiv:astro-ph/9910034. 
[13] M. Shibata and K. Taniguchi, Phys. Rev. D 73, 064027 

(2006), arXiv:astro-ph/0603145. 
[14] L. Baiotti, B. Giacomazzo, and L. Rezzolla, Phys. Rev. 

D78, 084033 (2008), 0804.0594. 
[15] L. Baiotti, B. Giacomazzo, and L. Rezzolla, Classical and 

Quantum Gravity 26, 114005 (2009), 0901.4955. 
[16] Y. T. Liu, S. L. Shapiro, Z. B. Etienne, and K. Taniguchi, 

Phys. Rev. D 78, 024012 (2008), 0803.4193. 
[17] B. Giacomazzo, L. Rezzolla, and L. Baiotti, MNRAS 399, 

L164 (2009), 0901.2722. 
[18] K. Kiuchi, Y. Sekiguchi, M. Shibata, and K. Taniguchi, 

Phys. Rev. D 80, 064037 (2009), 0904.4551. 
[19] M. D. Duez, Classical and Quantum Gravity 27, 114002 

(2010) , 0912.3529. 

[20] L. Rezzolla, L. Baiotti, B. Giacomazzo, D. Link, and J. A. 

Font, ArXiv e-prints (2010), 1001.3074. 
[21] M. Ruffert, M. Rampp, and H.-T. Janka, Astron. Astro- 

phys. 321, 991 (1997). 
[22] M. Ruffert, H.-T. Janka, and G. Schafer, Astron. Astro- 

phys. 311, 532 (1996). 
[23] M. Ruffert, H.-T. Janka, K. Takahashi, and G. Schaefer, 

Astron. Astrophys. 319, 122 (1997). 
[24] M. Ruffert and H.-T. Janka, Astron. Astrophys. 338, 535 

(1998). 

[25] R. Oechslin and H. Janka, MNRAS 368, 1489 (2006), 

arXiv:astro-ph/0507099. 
[26] M. Anderson, E. W. Hirschmann, L. Lehner, S. L. 

Liebling, P. M. Motl, D. Neilsen, C. Palenzuela, and J. E. 

Tohline, Phys. Rev. D 77, 024006 (2008), 0708.2720. 
[27] M. Anderson, E. W. Hirschmann, L. Lehner, S. L. 

Liebling, P. M. Motl, D. Neilsen, C. Palenzuela, and 

J. E. Tohline, Physical Review Letters 100, 191101 (2008), 

0801.4387. 

[28] S. Chawla, M. Anderson, M. Besselman, L. Lehner, S. L. 

Liebling, P. M. Motl, and D. Neilsen, Physical Review 

Letters 105, 111101 (2010), 1006.2839. 
[29] M. Ruffert and H. Janka, A&A 514, A66+ (2010). 
[30] M. Shibata and K. Uryu, Phys. Rev. D 74, 121503 (2006), 

arXiv:gr-qc/0612142. 



[31] M. Shibata and K. Uryu, Classical and Quantum Gravity 

24, 125 (2007), arXiv:astro-ph/0611522. 
[32] F. Foucart, M. D. Duez, L. E. Kidder, and S. A. Teukolsky, 

ArXiv e-prints (2010), 1007.4203. 
[33] M. Shibata, K. Kyutoku, T. Yamamoto, and K. Taniguchi, 

Phys. Rev. D 79, 044030 (2009), 0902.0416. 
[34] T. Yamamoto, M. Shibata, and K. Taniguchi, Phys. Rev. D 

78, 064054 (2008), 0806.4007. 
[35] Z. B. Etienne, Y. T. Liu, S. L. Shapiro, and T. W. Baum- 

garte, Phys. Rev. D 79, 044024 (2009), 0812.2245. 
[36] M. D. Duez, F. Foucart, L. E. Kidder, C. D. Ott, and S. A. 

Teukolsky, Classical and Quantum Gravity 27, 114106 

(2010), 0912.3528. 
[37] E. Rantsiou, S. Kobayashi, P. Laguna, and F. A. Rasio, 

ApJ 680, 1326 (2008), arXiv:astro-ph/0703599. 
[38] R. Popham, S. E. Woosley, and C. Fryer, ApJ 518, 356 

(1999), arXiv:astro-ph/9807028. 
[39] T. Piran, Reviews of Modern Physics 76, 1143 (2004), 

arXiv:astro-ph/0405503. 
[40] W. H. Lee and E. Ramirez-Ruiz, New Journal of Physics 

9, 17 (2007), arXiv:astro-ph/0701874. 
[41] P. Meszaros and M. J. Rees, MNRAS 257, 29P (1992). 
[42] N. Gehrels, E. Ramirez-Ruiz, and D. B. Fox, ARA&A 47, 

567 (2009), 0909.1531. 
[43] R. D. Blandford and R. L. Znajek, MNRAS 179, 433 

(1977). 

[44] V. V. Usov, Nature 357, 472 (1992). 
[45] B. D. Metzger, ArXiv e-prints (2010), 1001.5046. 
[46] J. Frank, A. King, and D. J. Raine, Accretion Power in As- 
trophysics: Third Edition (2002). 
[47] S. Rosswog and M. B. Davies, MNRAS 334, 481 (2002). 
[48] M. Lyutikov, ArXiv e-prints (2009), 0911.0349. 
[49] S. E. Woosley, A. Heger, and T. A. Weaver, Rev. Mod. 

Phys. 74, 1015 (2002). 
[50] Y. Sekiguchi and M. Shibata, Progress of Theoretical 

Physics 117, 1029 (2007), 0706.4154. 
[51] L. Dessart, C. D. Ott, A. Burrows, S. Rosswog, and 

E. Livne, ApJ 690, 1681 (2009), 0806.4380. 
[52] S. Rosswog and M. Liebendorfer, MNRAS 342, 673 

(2003), arXiv:astro-ph/0302301. 
[53] M. J. Rees and P. Meszaros, MNRAS 258, 41P (1992). 
[54] J. Goodman, A. Dar, and S. Nussinov, ApJ 314, L7 (1987). 
[55] M. Jaroszynski, Acta Astronomica 43, 183 (1993). 
[56] P. Meszaros, Rep. Prog. Phys. 69, 2259 (2006). 
[57] S. Kobayashi, B. Zhang, P. Meszaros, and D. Burrows, 

ApJ 655, 391 (2007), arXiv:astro-ph/0506157. 
[58] E. Nakar and T. Piran, ApJ 598, 400 (2003), arXiv:astro- 

ph/0303156. 

[59] B. Zhang, Y. Z. Fan, J. Dyks, S. Kobayashi, P. Meszaros, 
D. N. Burrows, J. A. Nousek, and N. Gehrels, ApJ 642, 
354 (2006), arXiv:astro-ph/0508321. 

[60] M. A. Abramowicz, M. Calvani, and L. Nobili, Nature 
302, 597 (1983). 

[61] J. C. B. Papaloizou and J. E. Pringle, MNRAS 208, 721 

(1984) . 

[62] J. C. B. Papaloizou and J. E. Pringle, MNRAS 213, 799 

(1985) . 

[63] Y. Kojima, Progress of Theoretical Physics 75, 251 (1986). 
[64] J. W. Woodward, J. E. Tohline, and I. Hachisu, ApJ 420, 
247 (1994). 

[65] J. A. Font and F. Daigne, Mon. Not. R. Astron. Soc. 334, 



32 



383 (2002). 

[66] O. Zanotti, L. Rezzolla, and J. A. Font, MNRAS 341, 832 

(2003), arXiv:gr-qc/0210018. 
[67] B. Paczynsky and P. J. Wiita, Astron. Astrophys. 88, 23 

(1980). 

[68] D. B. Wilson, Nature 312, 620 (1984). 

[69] M. A. Abramowicz, V. Karas, and A. Lanza, A&A 331, 

1143 (1998), arXiv:astro-ph/9712245. 
[70] F. Daigne and R. Mochkovitch, MNRAS 285, L15 (1997). 
[71] J. A. Font and F. Daigne, Astrophys.J 581, L23 (2002). 
[72] F. Daigne and J. A. Font, MNRAS 349, 841 (2004), 

arXiv:astro-ph/0311618. 
[73] R. Khanna and S. K. Chakrabarti, MNRAS 259, 1 (1992). 
[74] N. Masuda, S. Nishida, and Y. Eriguchi, MNRAS 297, 

1139 (1998). 

[75] S. Nishida, A. Lanza, Y. Eriguchi, and M. A. Abramow- 
icz, MNRAS 278, L41 (1996). 

[76] P. J. Montero, J. A. Font, and M. Shibata, Physical Review 
Letters 104, 191101 (2010), 1004.3102. 

[77] V. S. Safronov, Annales d'Astrophysique 23, 979 (1960). 

[78] A. Toomre, ApJ 139, 1217 (1964). 

[79] C. F. Gammie, ApJ 553, 174 (2001), arXiv:astro- 
ph/0101501. 

[80] Y. Kojima, Progress of Theoretical Physics 75, 1464 
(1986). 

[81] M. Abramowicz, M. Jaroszynski, and M. Sikora, A&A 
63, 221 (1978). 

[82] P. Goldreich, J. Goodman, and R. Narayan, MNRAS 221, 
339 (1986). 

[83] O. M. Blaes and W. Glatzel, MNRAS 220, 253 (1986). 
[84] R. Narayan, P. Goldreich, and J. Goodman, MNRAS 228, 
1 (1987). 

[85] J. Frank and J. A. Robertson, MNRAS 232, 1 (1988). 
[86] W. H. Zurek and W. Benz, ApJ 308, 123 (1986). 
[87] O. M. Blaes, MNRAS 227, 975 (1987). 
[88] O. M. Blaes and J. F. Hawley, ApJ 326, 277 (1988). 
[89] J. Goodman and R. Narayan, MNRAS 231, 97 (1988). 
[90] D. M. Christodoulou and R. Narayan, ApJ 388, 451 
(1992). 

[91] D. M. Christodoulou, ApJ 412, 696 (1993). 
[92] F. C. Adams, S. P. Ruden, and F. H. Shu, ApJ 347, 959 
(1989). 

[93] F. H. Shu, S. Tremaine, F. C. Adams, and S. P. Ruden, ApJ 
358, 495 (1990). 

[94] M. H. M. Heemskerk, J. C. Papaloizou, and G. J. 

Savonije, A&A 260, 161 (1992). 
[95] P. A. Taylor, J. C. Miller, and P. Podsiadlowski, ArXiv 

e-prints (2010), 1006.4624. 
[96] M. H. van Putten, Physical Review Letters 87, 091101 

(2001), arXiv:astro-ph/0107007. 
[97] T. Goodale, G. Allen, G. Lanfermann, J. Masso, T. Radke, 

E. Seidel, and J. Shalf, in Vector and Parallel Processing - 

VECPAR'2002, 5th International Conference, Lecture Notes 

in Computer Science (Springer, Berlin, 2003). 
[98] Cactus Website, Cactus Computational Toolkit, 

http: //www. cactuscode . org. 
[99] E. Schnetter, S. H. Hawley, and I. Hawke, Classical and 

Quantum Gravity 21, 1465 (2004), arXiv:gr-qc/0310042. 
[100] Carpet Website, Adaptive mesh refinement with Carpet, 

http: //www. carpetcode . org/. 
[101] E. Schnetter, P. Diener, E. N. Dorband, and M. Tiglio, 

Class. Quantum Grav. 23, S553 (2006), gr-qc/0602104. 
[102] B. Zink, E. Schnetter, and M. Tiglio, Phys. Rev. D 77, 

103015 (2008), 0712.0353. 



[103] E. Pazos, E. N. Dorband, A. Nagar, C. Palenzuela, 
E. Schnetter, and M. Tiglio, Classical and Quantum 
Gravity 24, 341 (2007), arXiv:gr-qc/0612149. 

[104] L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, and 
O. Rinne, Classical and Quantum Gravity 23, 447 (2006), 
arXiv:gr-qc/0512093. 

[105] J. Thornburg, Class. Quantum Grav. 4, 1119 (1987), URL 
|http: //stacks .io p . org/0264-9381 /4/1119| 

[106] J. Thornburg, Ph.D. thesis, University of British 
Columbia, Vancouver, British Columbia (1993). 

[107] R. Gomez, L. Lehner, R. Marsa, and J. Winicour, Phys. 
Rev. D 57, 4778 (1998), gr-qc/9710138. 

[108] R. Gomez, R. L. Marsa, and J. Winicour, Phys. Rev. D 56, 
6310 (1997), gr-qc/9708002. 

[109] R. Gomez, L. Lehner, R. Marsa, J. Winicour, A. M. Abra- 
hams, A. Anderson, P. Anninos, T. W. Baumgarte, N. T. 
Bishop, S. R. Brandt, et aL, Phys. Rev. Lett. 80, 3915 
(1998), gr-qc/9801069. 

[110] S. Bonazzola, E. Gourgoulhon, and J.-A. Marck, Pro- 
ceedings of the 19th Texas Symposium (1998), gr- 
qc/9904040. 

[Ill] L. E. Kidder, M. A. Scheel, and S. A. Teukolsky Phys. 

Rev. D 64, 064017 (2001), gr-qc/0105031. 
[112] E. Gourgoulhon, P. Grandclement, K. Taniguchi, 

J. Marck, and S. Bonazzola, Phys. Rev. D 63, 064029 

(2001). 

[113] P. Grandclement, S. Bonazzola, E. Gourgoulhon, and 
J.-A. Marck, J. Comput. Phys. 170, 231 (2001), gr- 
qc/0003072. 

[114] H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A. Teukol- 
sky, Computer Physics Communications 152, 253 (2003), 
arXiv:gr-qc/0202096. 

[115] J. Thornburg, Class. Quantum Grav. 21, 3665 (2004), gr- 
qc/0404059. 

[116] G. Calabrese and D. Neilsen, Phys. Rev. D 71, 124027 

(2005) , gr-qc/0412109. 

[117] M. A. Scheel, H. P. Pfeiffer, L. Lindblom, L. E. Kidder, 
O. Rinne, and S. A. Teukolsky, Phys. Rev. D 74, 104006 

(2006) , gr-qc/0607056. 

[118] H. P. Pfeiffer, D. Brown, L. E. Kidder, L. Lindblom, 
G. Lovelance, and M. A. Scheel, Class. Quantum Grav. 
24, S59 (2007), gr-qc/0702106. 

[119] F. Foucart, L. E. Kidder, H. P. Pfeiffer, and S. A. Teukol- 
sky, Phys. Rev. D 77, 124051 (2008), 0804.3787. 

[120] M. D. Duez, F. Foucart, L. E. Kidder, H. P. Pfeiffer, M. A. 
Scheel, and S. A. Teukolsky, Phys. Rev. D 78, 104015 

(2008) , 0809.0002. 

[121] E. Pazos, M. Tiglio, M. D. Duez, L. E. Kidder, and S. A. 

Teukolsky, Phys. Rev. D 80, 024027 (2009), 0904.0493. 
[122] P. C. Fragile, C. C. Lindner, P. Anninos, and J. D. 

Salmonson, ApJ 691, 482 (2009), 0809.3819. 
[123] M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. 

Matthews, and H. P. Pfeiffer, Phys. Rev. D 79, 024003 

(2009) , 0810.1767. 

[124] C. Ronchi, Journal of Computational Physics 124, 93 
(1996). 

[125] G. Calabrese and D. Neilsen, Phys. Rev. D 69, 044020 

(2004), gr-qc/0308008. 
[126] B. Szilagyi, D. Pollney, L. Rezzolla, J. Thornburg, and 

J. Winicour, Class. Quantum Grav. 24, S275 (2007), gr- 

qc/0612150. 

[127] C. Reisswig, N. T. Bishop, D. Pollney and B. Szilagyi, 
Classical and Quantum Gravity 27, 075014 (2010), 
0912.1285. 



33 



[128] L. D. Landau and E. M. Lifshitz, The Classical Theory of 

Fields, Course of Theoretical Physics, Volume 2 (Elsevier 

Butterworth-Heinemann, Oxford, 2004). 
[129] J. R. Wilson, Astrophys. J. 173, 431 (1972). 
[130] J. A. Font, Living Reviews in Relativity 6 (2003), URL 

http: //www. livingreviews . org/lrr- 2003-4 
[131] F. Banyuls, J. A. Font, J. M. Ibanez, J. M. Marti, and J. A. 

Miralles, Astrophys. J. 476, 221 (1997). 
[132] C. F. Gammie, J. C. McKinney and G. Toth, Astrophys. 

J. 589, 458 (2003), astro-ph/0301509. 
[133] P. Colella and P. Woodward, J. Comput. Phys. 54, 174 

(1984). 

[134] J. Marti and E. Mutter, J. Comput. Phys. 123, 1 (1996). 
[135] A. Harten, P. Lax, and B. van Leer, SIAM Rev. 25, 35 
(1983). 

[136] S. C. Noble, C. F. Gammie, J. C. McKinney, and L. Del 

Zanna, ApJ 641, 626 (2006), arXiv:astro-ph/0512420. 
[137] L. Lindblom and B. Szilagyi, Phys. Rev. D 80, 084019 

(2009), 0904.4873. 
[138] P. Secchi, Differential Integral Equations 9, 671 (1996), 

ISSN 0893-4983. 
[139] P. Secchi, Arch. Rational Me ch. Anal. 134, 155 (1996), 

ISSN 0003-9527, URL |http : //dx . doi . org/ 10 . 1007/| 

BF00379552I 

[140] J. Rauch, Trans. Amer. Math. Soc. 291, 167 (1985), ISSN 
0002-9947, URL h ttp : //dx . doi . org/10 . 2307/199 9902 

[141] T. P. Liu, "J. Differential Equations" 33, 92 (1979). 

[142] B. Gustafsson, H.-O. Kreiss, and J. Oliger, Time dependent 
problems and difference methods (Wiley New York, 1995). 

[143] L. Lehner, O. Reula, and M. Tiglio, Classical 
and Quantum Gravity 22, 5283 (2005), URL 
http: //www. citebase . org/abstract?id=oai : 

| arXiv.org:gr-qc/ 0507004 

[144] P. Diener, E. N. Dorband, E. Schnetter, and M. Tiglio, 
Journal of Scientific Computing 32, 109 (2007), URL doi : 
| 10.1007/sl0915-006-9123 r 7l 

[145] D. W. Neilsen and M. W. Choptuik, Classical and Quan- 
tum Gravity 17, 733 (2000), arXiv:gr-qc/9904052. 

[146] J. A. Font, H. Dimmelmeier, A. Gupta, and N. Ster- 
gioulas, Mon. Not. R. Astron. Soc. (2001), in press, astro- 
ph/0012477. 

[147] P. J. Montero, J. A. Font, and M. Shibata, ArXiv e-prints 

805 (2008), 0805.3099. 
[148] O. Korobkin, B. Aksoylu, M. Hoist, E. Pazos, and 

M. Tiglio, Class. Quantum Grav. 26, 145007 (2009), 

0801.1823. 

[149] B. Zink, O. Korobkin, E. Schnetter, and N. Stergioulas, 

Phys. Rev. D 81, 084055 (2010), 1003.0779. 
[150] N. Stergioulas and J. L. Friedman, ApJ 444, 306 (1995), 



arXiv:astro-ph/9411032. 
[151] S. Nishida and Y. Eriguchi, ApJ 427, 429 (1994). 
[152] H. Komatsu, Y. Eriguchi, and I. Hachisu, MNRAS 237, 

355 (1989). 

[153] M. Shibata, Phys. Rev. D 76, 064035 (2007). 

[154] J. D. Brown, Puncture evolution of schwarzschild black holes 

(2007), URL http : //www . citebase . org/abstract?id= 

oai: arXiv.org: 0705. 13 59 
[155] S. W. Andalib, J. E. Tohline, and D. M. Christodoulou, 

ApJS 108, 471 (1997). 
[156] D. Brown, P. Diener, O. Sarbach, E. Schnetter, and 

M. Tiglio, Phys. Rev. D79, 044023 (2009), 0809.3533. 
[157] H. A. Williams and J. E. Tohline, ApJ 315, 594 (1987). 
[158] J. A. Font and F. Daigne, MNRAS 334, 383 (2002), 

arXiv:astro-ph/0203403. 
[159] C. W. Misner, K. S. Thorne, and J. A. Wheeler, Gravitation 

(W. H. Freeman, San Francisco, 1973). 
[160] E. Schnetter, Class. Quantum Grav. 20, 4719 (2003), gr- 

qc/0306006, URL http : //stacks . iop . org/0264-9381/ 

20/47191 

[161] A. T. Okazaki, S. Kato, and J. Fukue, PASJ 39, 457 (1987). 
[162] L. Rezzolla, S. Yoshida, and O. Zanotti, MNRAS 344, 978 

(2003), arXiv:astro-ph/0307488. 
[163] H. Dimmelmeier, N. Stergioulas, and J. A. Font, MNRAS 

368, 1609 (2006), arXiv:astro-ph/0511394. 
[164] J. A. Font, T. Goodale, S. Iyer, M. Miller, L. Rezzolla, 

E. Seidel, N. Stergioulas, W.-M. Suen, and M. Tobias, 

Phys. Rev. D 65, 084024 (2002), gr-qc/0110047. 
[165] L. Baiotti, I. Hawke, P. J. Montero, F. Loffler, L. Rezzolla, 

N. Stergioulas, J. A. Font, and E. Seidel, Phys. Rev. D 71, 

024035 (2005), gr-qc/ 0403029. 
[166] E. B. Abdikamalov, H. Dimmelmeier, L. Rezzolla, and 

J. C. Miller, MNRAS 392, 52 (2009), 0806.1700. 
[167] C. D. Ott, Classical and Quantum Gravity 26, 063001 

(2009), 0809.0695. 
[168] E. B. Abdikamalov, C. D. Ott, L. Rezzolla, L. Dessart, 

H. Dimmelmeier, A. Marek, and H. Janka, Phys. Rev. D 

81, 044012 (2010), 0910.2703. 
[169] S. Brandt and E. Seidel, Phys. Rev. D 54, 1403 (1996). 
[170] R. P. Kerr and A. Schild, General Relativity and Gravita- 
tion 41, 2485 (2009). 
[171] R. Takahashi, MNRAS 382, 567 (2007). 
[172] G. Calabrese, L. Lehner, O. Reula, O. Sarbach, and 

M. Tiglio, Class. Quantum Grav. 21, 5735 (2004), gr- 

qc/0308007. 

[173] J. A. Font, N. Stergioulas, and K. D. Kokkotas, Mon. Not. 
R. Astron. Soc. 313, 678 (2000). 



